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Point  defects,  submicroscopic  clusters  and  extended  defects  play  an  important 
role  in  dopant  diffusion  in  ion-implanted  silicon.  Point  defects  aggregate  during  the 
annealing  process  to  form  small  clusters,  {311}  defects  and  dislocation  loops.  {311} 
defects,  in  turn,  have  been  observed  to  release  interstitials  upon  annealing  the  silicon. 
These  interstitials  interact  with  dopant  atoms  such  as  boron  leading  to  transient  enhanced 
diffusion. 

In  this  work,  molecular  dynamics  simulations  are  used  to  study  the  motion  of 
interstitial,  vacancy  and  di-interstitial  species  in  silicon.  From  these  calculations  the 
validity  of  the  new  environment  dependent  interatomic  potential  is  determined  by 
comparing  the  interstitial  migration  energy  barrier  obtained  from  different  interatomic 
potentials. 

Molecular  dynamics  simulations  are  also  used  to  determine  the  formation  energies 
of  self-interstitial  clusters  with  less  than  10  interstitials.  By  comparing  these  formation 
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energies  with  those  of  { 3 1 1 } defects  of  the  same  size  inferences  are  drawn  about  the 
transition  from  these  sub-microscopic  clusters  to  {31 1 } defects.  Further,  in  this  work,  the 
formation  energies  of  {31 1 } defects  and  perfect  dislocation  loops  are  computed  using  the 
conjugate  gradient  technique.  From  the  formation  energy  values  of  {311}  defects  with 
varying  widths,  the  possibility  of  a length  to  width  correlation  of  these  defects  is  studied 
which  has  a bearing  on  the  traditional  view  of  the  growth  of  these  defects. 

Finally,  the  relative  stability  of  {31 1}  defects  and  perfect  dislocation  loops  is 
determined  from  the  nucleation  free  energies  of  these  defects.  From  this  study  a better 
understanding  of  the  unfaulting  behavior  of  {31 1}  defects  to  form  dislocation  loops  is 
gained. 
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CHAPTER  1 
INTRODUCTION 

1.1  Motivation 

In  the  semiconductor  industry,  ion  implantation  is  used  most  commonly  for 
producing  reproducible,  precisely  controlled  doping  in  the  fabrication  of  integrated 
circuits.  Ion  implantation  is  also  used  to  modify  the  properties  of  solids.  This  work  will 
deal  with  the  application  of  ion  implantation  in  semiconductors,  primarily,  for  doping  of 
silicon.  The  major  disadvantage  of  ion  implantation  is  that  the  implant  produces  damage 
in  silicon.  This  necessitates  a post-implant  annealing  step  to  activate  the  dopants.  As  a 
result  of  this  anneal,  silicon  self-interstitials  formed  by  the  implant  attain  sufficient 
energy  to  aggregate  into  extended  defects.  These  defects,  by  acting  as  a source  of 
interstitials,  lead  to  anomalous  diffusion  of  dopants  due  to  the  interaction  between  self- 
interstitials and  dopants  upon  elevated  temperature  annealing.  This  phenomenon,  known 
as  transient-  enhanced  diffusion  (TED),  results  in  changes  in  the  junction  depth  leading  to 
degradation  in  the  performance  of  future  generation  semiconductor  devices. 

TED  also  is  found  to  occur  under  conditions  where  there  is  a suppression  of  the 
formation  of  extended  defects  suggesting  that  submicroscopic  interstitial  defects  exist 
[Zhang  95].  While  structural  characterization  of  {311}  defects  is  possible  using 
Transmission  electron  microscopy  (TEM),  the  characteristics  of  small  clusters  are  not 
well  understood.  It  was  only  recently  that  evidence  of  the  presence  of  small  clusters  was 
found  for  self-implants  using  deep  level  transition  spectroscopy  [Benton  97],  Moreover, 
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for  many  years  now,  considerable  experimental  and  theoretical  investigations  have  been 
carried  out  to  study  the  structure,  energetics  and  migration  of  point  defects  in  silicon  in 
addition  to  their  role  in  TED.  Experimental  work  in  silicon  has  produced  conflicting 
evidence  resulting  in  debates,  particularly,  on  the  relative  contribution  of  vacancies  and 
interstitials  on  silicon  self-diffusion. 

In  contrast  to  this,  theoretical  work  has  made  significant  progress  in  this  area 
using  both  first  principle  calculations  and  atomistic  simulations  based  on  empirical 
interatomic  potentials.  However,  controversy  remains  regarding  some  of  the  most 
fundamental  properties  of  point  defects.  A number  of  crucial  parameters  such  as  the 
diffusivities  and  formation  energies  of  point  defects  have  not  been  well  established. 

Moreover,  small  interstitial  clusters  have  not  been  sufficiently  characterized  and 
the  transition  from  small  clusters  to  extended  defects  needs  further  investigation.  On  the 
issue  of  extended  defects  such  as  {311}  defects  and  dislocation  loops,  parameters  such  as 
formation  energies  and  relative  stability  of  these  extended  defects  have  not  been 
conclusively  obtained.  Some  work  has  been  done  related  to  Si  self-diffusion  [Gilmer  95] 
using  the  Stillinger-Weber  interatomic  potential.  As  far  as  small  silicon  self-interstitial 
clusters  are  concerned,  tight  binding  molecular  dynamics  calculations  have  been  done  to 
obtain  information  about  mainly  the  structure  of  clusters  with  a maximum  of  nine 
interstitials  [Colombo  99],  More  recently,  the  transition  from  small  interstitial  clusters  to 
extended  {311}  defects  has  been  investigated  experimentally  using  photoluminescence 
and  TEM  studies  [Coffa  2000],  The  conclusion  of  this  study  was  that  small  interstitial 
clusters  are  not  direct  precursors  of  the  {311}  defects  and  that  a structural  transformation 
must  occur  at  some  stage  of  the  growth  process. 
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Nucleation  of  vacancy  type  dislocation  loops  and  voids  has  been  studied  [Tan  97], 
However  there  is  no  known  work  in  the  literature  regarding  the  interstitial-type 
dislocation  loops.  Moreover,  the  relative  stability  of  loops  versus  (31 1 }s  has  not  been 
studied.  The  unfaulting  of  {31 1}  defects  has  been  found  to  be  the  source  of  dislocation 
loops  by  TEM  studies  [Li  98].  Understanding  why  {3 1 1 }s  form  rods  is  essential  to 
developing  accurate  {311}  nucleation,  growth,  dissolution  models  and  loop  nucleation 
models. 


1 .2  Point  Defect  Diffusion  in  Silicon 

With  shrinking  device  dimensions,  it  becomes  important  to  understand  the 
tendency  of  dopants  to  exhibit  enhanced  diffusion.  Vacancies  and  interstitials  produced 
due  to  implant  damage  have  a deleterious  effect  on  the  junction  depth  produced  during 
post-implant  annealing.  As  an  example,  release  of  the  interstitials  by  {311}  defects  tends 
to  produce  anomalous  enhanced  diffusion  by  the  interaction  of  these  interstitials  and 
dopant  atoms  such  as  boron  or  phosphorus.  Excess  interstitials  are  generated  by  a “knock 
on”  effect.  Implanted  ions  displace  atoms  from  their  lattice  sites  close  to  the  surface  and 
create  a shallow  vacancy-rich  region.  The  displaced  atoms  and  the  implanted  ions  come 
to  rest  deeper  inside  the  silicon  creating  an  interstitial-rich  region  there. 

At  present,  it  is  believed  that  most  of  the  Frenkel  pairs  recombine  and  leave  an 
interstitial  concentration  approximately  equal  to  the  dose.  The  aggregation  of  these 
excess  interstitials  leads  to  the  formation  of  {31 1 } defects  and  loops  which  in  turn  act  as 
a source  of  interstitial  emission  and  thereby,  enhanced  diffusion. 

Large  leakage  currents  are  seen  near  p-n  junction  depletion  regions  due  to  the 
presence  of  defects  in  that  area.  Moreover,  there  is  a competition  between  self-interstitials 
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and  inactive  dopants  for  substitutional  lattice  sites.  This  affects  the  electrical  activation  of 
implanted  dopants. 

It  has  been  observed  that  there  are  many  mechanisms  for  the  diffusion  of  point 
defects  through  a crystal  lattice.  The  most  common  mechanisms  are  described  below. 

• Vacancy  mechanism:  This  occurs  when  a dopant  atom  moves  to  occupy  a lattice  site 
and  thus  becomes  a substitutional  atom.  In  order  to  achieve  long-range  movement,  the 
vacancy,  which  was  originally  in  the  dopant  site  must  diffuse  to  a location  some 
distance  away  from  the  original  lattice  site.  A vacancy  must  diffuse  at  least  three 
lattice  sites  away  in  order  to  complete  the  diffusion  step.  This  happens  to  be  the 
minimum  distance  that  allows  the  vacancy  to  come  back  to  a site  adjacent  to  the 
dopant  atom  through  another  path  [Fahey  89]. 

• Interstitialcy  diffusion:  This  occurs  when  a silicon  interstitialcy  pair  approaches  a 
substitutional  dopant  atom’s  lattice  site.  Diffusion  occurs  when  one  of  the 
interstitialcy  -silicon  atom  forms  a bond  with  a third  silicon  atom  to  produce  a new 
interstitialcy  pair.  A similar  mechanism  occurs  with  a silicon-dopant  interstitialcy 
defect.  It  should  be  noted  that,  unlike  the  vacancy  mechanism  here  the  diffusing 
defect  remains  intact  [Fahey  89], 

• Interstitial  mechanism:  This  mechanism  is  operative  when  an  interstitial  silicon  atom 
goes  near  a substitutional  dopant  atom,  dislodges  it  from  its  site  and  occupies  the 
lattice  site.  Therefore  this  mechanism  is  called  a “kickout  mechanism”.  The  impurity 
atom  is  now  liberated  and  can  freely  migrate  in  the  lattice  interstices,  until  it  takes 
part  in  another  interstitial  exchange  [Fair  81,  Fahey  85, Tan  85],  Fig.  1.1  illustrates 
these  three  types  of  dopant  diffusion  mechanisms. 
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1 .3  Transient  Enhanced  Diffusion 

For  many  years  to  come,  ion  implantation  will  continue  to  be  the  mainstay  of 
introducing  dopants  into  silicon  during  the  fabrication  of  devices.  During  implantation, 
damage  is  produced  in  the  form  of  vacancies  and  interstitials.  In  order  to  remove  the 
damage  the  implanted  silicon  wafer  is  subjected  to  a high  temperature  anneal.  However, 
during  this  annealing  phase,  many  types  of  extended  defects  are  formed.  During  this 
anneal  anomalous  diffusion  of  dopants  is  found  to  occur  which  is  different  from 
equilibrium  diffusion  as  postulated  by  Fick’s  law. 

Specifically,  the  diffusion  of  boron  is  enhanced  several  orders  of  magnitude 
greater  than  the  diffusion  of  boron  in  thermal  equilibrium  [Michel  87],  The  TED  of  boron 
takes  place  primarily  because  of  silicon  excess  self-interstitials  which  are  created  during 
boron  implantation.  From  theoretical  calculations,  such  as  first  principle  calculations,  the 
activation  energy  of  boron  diffusion  associated  with  Si  interstitials  is  the  lowest  among 
possible  diffusion  mechanisms  [Car  84,  Nichols  89],  This  enhanced  diffusion  is  transient 
in  that  it  occurs  for  a certain  limited  time  period,  which  depends  on  the  annealing 
temperature.  Hence,  this  anomalous  diffusion,  which  is  many  orders  of  magnitude  larger 
than  equilibrium  diffusion,  is  called  transient  enhanced  diffusion  (TED). 

TED  is  found  to  be  the  result  of  a complicated  interaction  between  dopant  atoms, 
point  defects  and  extended  defects.  Investigations  [Eaglesham  94]  have  shown  that 
extended  defects  such  as  the  {311}  defects  can  act  as  sources  of  interstitials  for  driving 
TED.  It  has  been  suggested  that  Si  interstitials,  which  are  released  from  extended  defects, 
aid  dopant  impurity  movement  via  the  “kickout  mechanism”  discussed  earlier.  Other 
investigations  [Zhang  95]  have  suggested  that  TED  can  occur  even  in  the  absence  of 
{311}  formation  due  to  Boron  interstitial  clusters,  which  act  as  a source  of  interstitials. 
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1 .4  Classification  of  Extended  Defects  in  Silicon 

During  the  post-implant  annealing  of  the  silicon  different  types  of  extended 
defects  form.  The  type  of  extended  defect  formed  depends  on  the  implant  and  anneal 
conditions.  A systematic  analysis  of  these  defects  has  been  done  leading  to  a 
classification  scheme  [Jones  87],  A schematic  of  the  connection  between  the  damage 
density  distribution  and  the  effective  threshold  damage  density  is  illustrated  in  fig.  1.2. 

1.4.1  Type  I:  Sub-Amorphization  Defects 

This  damage  is  observed  when  the  dose  exceeds  a critical  value  and  there  is  no 
concomitant  amorphization.  These  defects  form  around  the  projected  range  of  the 
implanted  ions.  These  defects  have  two  types  of  morphology:  rod-like  {311}  defects  and 
dislocation  loops.  {3 1 1 }s  form  at  lower  doses  than  dislocation  loops.  The  threshold  dose 
has  been  shown  to  be  as  low  as  7xl012/cm2  whereas  the  threshold  for  dislocation 
formation  is  around  2x1 014/  cm2.  These  subamorphization  defects  form  due  to  the  fact 
that  the  concentration  of  ions  implanted  exceeds  the  equilibrium  vacancy  concentration. 
Upon  Frenkel  pair  recombination  the  implantation  process  results  in  a supersaturation  of 
interstitials.  When  the  implant  dose  exceeds  2x1 014/  cm2  but  below  the  amorphization 
threshold  then  both  {311}  defects  and  dislocation  loops  form.  If  the  anneal  is  above  800 
°C  after  implanting  with  a dose  below  2xl014/  cm2  then  no  extended  defects  are  seen  after 
annealing.  For  higher  doses  without  amorphization  then  subthreshold  dislocation  loops 
are  observed  which  can  be  stable  to  temperatures  above  1000  °C. 

1.4.2  Type  II:  End-of-Range  Defects 

These  defects  arise  when  an  amorphous  layer  is  produced.  The  location  of  these 
defects  is  below  the  amorphous-crystalline  interface  in  heavily  damaged  but  crystalline 
material.  These  defects  consist  of  both  dislocation  loops  and  {311}  defects  at  low 
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annealing  temperature.  However,  unlike  the  subthreshold  defects,  where  there  is  a 
possibility  of  formation  of  {31 1 } defects  without  loops,  the  concentration  of  interstitials 
beyond  the  amorphous/crystalline  interface  is  such  that  both  defects  are  observed 
together.  There  are  two  major  sources  of  the  interstitials  for  these  defects. 

The  first  source  consists  of  transmitted  ions  which  rest  below  the 
amorphous/crystalline  interface.  The  second  source  is  the  recoiling  of  the  excess 
interstitials  deeper  into  the  material  because  of  the  forward  momentum  of  the  ion  beam. 
The  shallower  region  where  many  of  the  excess  vacancies  reside  in  non-amorphizing 
implants  is  now  amorphized.  Therefore,  after  solid  phase  regrowth  there  are  excess 
interstitials  from  the  separated  Frenkel  pairs  which  could  not  recombine  because  of 
vacancy  incorporation  during  regrowth.  These  excess  interstitials  can  now  constitute  the 
end  of  the  end-of  range  loop  density  as  well  as  providing  the  general  interstitial 
supersaturation  in  the  end-of  range  region.  The  concentration  of  type  II  defects  is  not  as 
strongly  dependent  on  the  dose  as  in  the  case  of  type  I defects.  This  is  because  an 
amorphous  layer  exists  which  increases  in  thickness  as  dose  increases. 

High  annealing  temperatures  are  necessary  to  dissolve  end-of  range  defects, 
particularly,  the  loops  completely.  These  temperatures  are  not  attained  in  the  thermal 
budgets  of  current  processing  techniques.  If  these  loops  exist  in  the  space  charge  regions 
of  junctions  they  can  cause  high  leakage  currents.  Therefore,  the  junction  should  be 
located  deep  enough  such  that  the  end-of  range  defects  are  not  in  the  space  charge  region 
of  the  device. 

1.4.3  Type  III  Defects 

Formation  of  this  category  of  defects  can  occur  due  to  regrowth  of  the  amorphous 
layer.  These  defects  are  found  to  exist  in  the  region  previously  occupied  by  the 
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amorphous  layer.  There  are  many  forms  of  this  defect  type.  Hairpin  dislocations  form 
when  the  regrowing  amorphous/crystalline  interface  encounters  pockets  of  misoriented 
crystalline  material.  The  other  form  of  type  III  defects  are  microtwins.  Defect  formation 
also  occurs  just  beneath  the  surface.  This  is  due  to  the  rejection  of  some  species,  such  as 
flourine,  by  the  progressing  amorphous/crystalline  interface  for  sufficient  dose. 

1 .4.4  Type  IV  Defects 

These  defects  form  due  to  the  fact  that  regrowth  of  a buried  amorphous  layer 
occurs  by  the  movement  of  both  amorphous/crystalline  interfaces.  At  the  depth  where 
these  two  interfaces  meet  type  IV  defects  originate.  These  defects  are  also  called  “zipper” 
or  “clamshell”  defects. 

1 .4.5  Type  V Defects 

In  ion  implantation,  which  is  a non-equilibrium  process,  the  solid  solubility  of  a 
species  can  be  exceeded  in  a target  material.  During  solid  phase  epitaxy,  dopant 
concentrations  in  excess  of  the  solid  solubility  can  be  incorporated  in  the  lattice  sites  as 
substitutional  species.  Upon  further  annealing,  precipitation  occurs.  The  defects  thus 
formed  at  a dopant  projected  range  depth  consist  of  both  dislocation  loops  and 
precipitates. 

1 .5  Extended  1311)  Defects 

The  {311}  extended  defect  is  peculiar  to  silicon  and  germanium.  Mathews  and 
Ashby  [Mathews  73]  observed  these  defects  in  irradiated  n-type  and  p-type  silicon  at 
temperatures  above  400  °C.  They  found  these  defects  to  be  similar  in  electron-irradiated 
and  proton-irradiated  silicon.  {311}  defects  were  for  the  first  time  found  in  electron- 
irradiated  germanium  by  Ferreira  et  al.  [Ferreira  76]  , and  in  annealed,  ion-implanted 
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silicon  by  Davidson  et  al.  and  Chadderton  et  al.  [Davidson  70,  Chadderton  71],  It  is  a 
{311}  stacking  fault  since  fringe  contrast  is  similar  to  that  observed  in  the  common 
{111}  stacking  fault.  However,  the  defect  has  {311}  habit  planes.  Initially,  the  only 
structural  model  proposed  was  that  there  is  a layer  of  interstitial  atoms  in  these  defects 
each  with  three  dangling  bonds  with  bonds  of  the  matrix  crystal  atoms  on  either  side  of 
the  interstitial  atom  layer  stretched  by  35-48%  [Salisbury  79],  This  structure  is  not 
favorable  energetically. 

Tan  [Tan  81]  proposed  a low-energy  model  for  {311}  defects  which  consists  of  a 
layer  of  interstitial  atoms  with  no  dangling  bonds  and  with  insignificant  bond  length 
changes.  This  model  is  shown  in  fig.  1 .3  (a)  and  fig.  1 .3  (b).  In  fig.  1 .3  (a)  the  interstitial 
atoms  are  on  the  (1  1 0)  plane  layers  of  the  matrix  crystal.  The  bond  length  changes  are 
within  5%  for  all  atoms.  Figure  1 .3  (b)  shows  that  excepting  for  the  single  bond  which  is 
shared  by  two  neighboring  seven-membered  rings  and  is  stretched  by  1 5%,  all  the 
remaining  bonds  are  more  or  less  unchanged  in  their  lengths.  In  both  these  depictions  the 
defect  habit  plane  is  ( 1 1 3). 

{311}  defects  are  a part  of  subthreshold  or  type  I defects.  As  mentioned  in  a 
previous  section,  the  dissolution  of  these  {311}  defects  is  one  of  the  principal  sources  of 
interstitials  causing  TED.  Below  the  loop  formation  threshold  the  interstitials  generated 
due  to  the  implantation  are  in  a supersaturated  state  upon  Frenkel  pair  recombination  and 
these  interstitials  agglomerate  to  form  the  {311}  defects.  The  {311}  defects  are  unstable 
compared  to  the  loops  and  they  dissolve  upon  annealing  the  silicon  between  700  and  800° 
C.  Figure  1 .4  shows  the  loss  of  interstitials  from  {311}  defects  as  a function  of  annealing 
time  at  750°  C [Jones  96], 
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As  mentioned  in  the  previous  section,  type  II  defects  or  end  of  range  defects  form 
when  the  silicon  has  been  amorphized  during  implantation.  Both  {311}  defects  and 
dislocation  loops  are  seen  to  form  just  below  the  amorphous/crystalline  interface  at  low 
annealing  temperatures.  Unlike  in  subthreshold  defects,  where  {311}  defects  can  form 
without  loops,  in  the  end  of  range  defects  both  loops  and  {3 1 1 }s  are  seen  to  form.  The 
{311}  defects  which  are  present  in  the  end  of  range  regime  coarsen  and  dissolve  upon 
annealing  the  sample  between  700  and  800°  C.  During  the  dissolution  of  the  {31 1 }s  TED 
is  seen  to  occur  for  dopants  located  both  below  the  amorphous  crystalline  interface  and  in 
the  regrown  region. 

The  atomic  structure  and  energy  of  {31 1 } defects  have  been  studied  by  energy 
minimization  calculations  [Kohyama  92]  although  the  effect  of  the  width  of  the  defect  on 
the  energetics  of  the  defect  has  not  been  completely  investigated.  Using  the  Tersoff 
interatomic  potential  the  energetics  of  {31 1 } defects  has  been  investigated,  for  infinitely 
long  chains,  as  a function  of  defect  size  and  width  [Cuendet  96].  Continuum  calculations 
have  been  performed  in  the  modeling  of  {31 1 } defects  [Hobler  99],  Also,  total  energy 
calculations  using  the  tight-binding  scheme  have  been  carried  out  to  study  the  structural 
properties  and  energetics  of  {31 1}  defects  [Kim  97]. 

1 .6  Scope  of  This  Study 
The  scope  of  this  work  is  as  follows: 

• To  carry  out  a comparative  study  of  diffusivities  of  point  defects  using  Molecular 
Dynamics  simulations  with  three  different  interatomic  potentials,  namely,  the 
Stillinger  Weber  (SW)  potential,  the  Tersoff  potential  and  the  recently  developed 
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Environment  Dependent  Interatomic  potential  (EDIP).  From  this  study  the  validity  of 
the  EDIP  potential  is  determined. 

• To  calculate  formation  energies  of  small  interstitial  clusters,  using  Molecular 
Dynamics  simulations  with  the  SW  potential,  and  to  analyze  the  transformation  of 
these  clusters  into  {311}  defects. 

• To  calculate  the  formation  energies  of  {311}  defects,  using  the  conjugate  gradient 
technique,  as  a function  of  size  and  width  and  to  calculate  the  formation  energy  of 
perfect  dislocation  loops,  also  using  the  conjugate  gradient  method.  Thereby, 
inferences  are  drawn  about  the  relative  stability  of  these  two  types  of  extended 
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Fig  1.1  Dopant  diffusion  mechanisms 
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Fig.  1 .2  Schematic  of  damage  classification  after  ion  implantation  and  subsequent 
annealing 
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(a)  Atomic  model  of  a (3 1 1 } defect  with  a matrix  crystal  displacement  of 

(b)  Atomic  model  of  a {311}  defect  with  a matrix  crystal  displacement  of 


1/5 

1/25 


Trapped  Interstitials  in  {311}'s  (cm 


15 


Fig.  1.4  Interstitial  loss  from  {311}  defects  as  a function  of  annealing  time  at  750°  C 
obtained  using  quantitative  PTEM.  Implant  conditions  are  20  KeV  B+,  2 x 1014  /cm2 


CHAPTER  2 

THEORY  OF  COLLISIONS 
2.1  Ion  Implantation  Models 

Ion  implantation  models  can  be  either  analytical  theories  (AT)  or  computer 
simulations  (CS).  According  to  both  the  models  the  implanted  ion  loses  energy  by  elastic 
and  inelastic  collisions  [Townsend  76].  When  ions  collide  with  the  atomic  nuclei  the 
collision  is  elastic  and  collisions  between  ions  and  electrons  are  inelastic  [Bohr  48],  The 
losses  due  to  inelastic  collisions  are  simple  energy  losses  [Lindhard  61].  The  collision 
will  not,  in  this  case,  lead  to  a deviation  of  the  trajectory.  On  the  other  hand,  elastic 
collisions  modify  the  ion’s  trajectory  and  the  scattering  angle  in  the  collision  is 
determined  by  an  interatomic  potential.  The  binary  collision  approximation,  which  is  the 
simplest  model  of  collisions,  assumes  that  the  collision  takes  place  only  between  two 
atoms  in  the  solid.  The  distance  between  colliding  atoms  alone  is  used  in  implementing 
the  interatomic  potential.  Both  analytical  theories  and  computer  simulations  employ  such 
interatomic  potentials.  As  a result,  these  models  describe  the  damage  profiles  and  range 
of  implanted  ions  in  solids  as  well  as  the  sputtering  of  atoms  from  the  surface  of  the  solid. 

Since  these  potentials  do  not  account  for  the  way  atoms  are  bound  in  the  solid, 
more  advanced  potentials  are  used  to  determine  the  interatomic  interactions  in  a solid. 
Common  examples  of  these  potentials  are  the  EAM  [Daw  84]  potential  for  metals  and  the 
Tersoff  [Tersoff  89]  and  Stillinger  and  Weber  (SW)  [Stillinger  85]  potentials  for  covalent 
materials  like  silicon.  The  SW  potential  consists  of  two  and  three  body  terms  and  was 
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fitted  to  experimental  properties  of  the  diamond  cubic  and  molten  phases  of  silicon. 
Among  its  uses  are  the  study  of  lattice  dynamics  [Li  88],  point  defects  [Batra  87], 
surfaces  [Poon  90],  and  the  liquid  and  amorphous  phases  [Stillinger  85,  Broughton  87], 
[Luedtke  88,  Kluge  87,  Ishimaru  96],  The  Tersoff  potential  contains  many-body 
interactions  included  in  a bond  order  term  fitted  to  various  silicon  polytypes.  Its  uses 
include  the  study  of  lattice  dynamics  [Li  88],  point  defects  [Tersoff  88], 
thermomechanical  properties  [Porter  97]  and  the  liquid  and  amorphous  phases  [Tersoff 
88,  Ishimaru  96,  Kelires  88], 

Recently,  more  than  30  empirical  interatomic  potentials  have  been  developed  for 
silicon.  These  potentials  differ  in  the  degree  of  sophistication,  functional  form  , and  range 
of  interaction.  No  single  potential  can  accurately  model  various  special  atomic 
configurations.  It  has  been  claimed  that  the  transferability  of  the  SW  and  Tersoff 
potentials  to  a wide  class  of  structures  remains  in  doubt  [Justo  98],  It  has  also  been 
claimed  that,  the  most  recent  interatomic  potential  called  the  environment-dependent 
interatomic  potential,  provides  a better  description  of  local  structures  and  is,  at  the  same 
time,  computationally  efficient. 

These  interatomic  potentials  can  be  used  in  molecular  dynamics  to  determine  the 
evolution  of  a group  of  atoms  as  a function  of  time.  This  is  done  by  solving  Newton’s 
equations  of  motion  involving  the  potentials  to  calculate  energies  and  forces  of  each  atom 
in  the  computational  cell.  Although  multiple  collisions  between  target  atoms  are  included 
in  the  calculations,  electronic  effects  are  not  included.  Certain  models  like  tight  binding 
combine  the  advantages  of  classical  molecular  dynamics  (using  interatomic  potentials) 
and  ab  initio  calculations.  Classical  molecular  dynamics  is,  by  far,  the  best  technique  for 
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studying  ion  implantation  in  an  energy  range  between  eV  and  several  keV.  By  employing 
a sound  interatomic  potential  in  classical  molecular  dynamics  calculations  it  would  be 
possible  to  gain  insight  into  trends  and  processes  involved  in  ion  implantation  collisions, 
something  which  is  impossible  using  analytical  models  and  binary  collision  simulations. 

2.2  Ion  Range  and  Stopping  in  Solid  Materials 

The  observation  of  emission  of  energetic  particles  from  radioactive  materials 
created  interest  in  how  these  particles  were  slowed  down  in  matter.  In  the  early  days, 
since  there  was  no  accurate  model  for  the  atom,  attempts  to  create  a particle  energy  loss 
theory  were  futile.  J.J. Thomson  was  the  first  to  address  the  issue  of  scattering  of  two 
point  charges  in  his  classic  book  on  electricity  [Thomson  03],  Although  this  book 
contains  an  excellent  description  of  classical  Coulombic  scattering  between  energetic 
charged  particles,  it  does  not  calculate  actual  stopping  powers. 

In  the  decade  following  the  discovery  of  radioactivity,  sufficient  experimental 
evidence  was  collected  to  make  stopping  power  theory  one  of  the  dominant  issues  for 
those  trying  to  develop  a model  for  the  atom.  Both  J.J.Thomson  and  Niels  Bohr  were 
actively  involved  in  the  analysis  of  the  stopping  of  charged  particles  by  matter.  However, 
the  model  of  the  atom  as  an  entity  with  a heavy  positively-charged  core  is  attributed  to 
Bohr’s  work  [Bohr  13,  Bohr  15].  Bohr’s  initial  work  is  informative  because,  for  the  first 
time,  a unified  theory  of  stopping  was  proposed.  One  of  Bohr’s  original  postulates  was 
that  the  energy  loss  of  ions  passing  through  matter  has  two  contributors  namely,  nuclear 
stopping  (energy  loss  to  the  material’s  positive  atomic  cores)  and  electronic  stopping 
(energy  loss  to  the  material’s  light  electrons). 
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The  Lindhard,  Sharff  and  Schiott  theory  (LSS)  developed  in  1963  was  the  first 
analytical  approximation  to  the  range  and  stopping  of  ions  [Lindhard  63],  According  to 
the  terminology  introduced  by  these  authors,  a reduced  energy  e and  a reduced  length  p 
as  follows: 

aM, 

e = E — (2  i t 

Z^e2^  + m2)  K ' ' 
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47ta2M1 
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Subscript  1 denotes  the  projectile  and  subscript  2 refers  to  the  target  atom,  M 
refers  to  the  mass  and  Z the  atomic  number.  The  existence  of  electrons  in  an  atom  is 
accounted  for  in  the  screening  distance  a.  Q is  the  target  density.  The  energy  loss  of  an 
ion,  which  is  interacting  with  the  target,  can  be  equated  to  the  sum  of  elastic  or  nuclear 
collisions  and  the  inelastic  or  electronic  collision.  This  equation  can  be  represented  as 
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where  the  subscript  n refers  to  nuclear  and  subscript  e refers  to  electronic.  sn  and  se  are 
called  elastic  and  inelastic  stopping  powers  and  these  can  be  obtained  analytically 
assuming  very  simple  interatomic  potentials  and  random  collisions.  Figure  2.1  shows  the 
nuclear  and  electronic  stopping  powers  as  a function  of  reduced  energy  obtained  using 
the  assumptions  of  the  LSS  theory  [Robinson  94], 

Whereas  the  elastic  energy  loss  is  independent  of  the  mass  ratios,  the  inelastic 


energy  loss  is  more  important  for  high  mass  ratios.  Upon  comparing  the  contribution  of 
the  two  types  of  collisions,  it  is  seen  that  for  low  values  of  £ the  elastic  collisions 
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dominate  and  for  e>l  the  inelastic  collisions  become  significant.  It  should  also  be  noted 
that  even  at  low  energies,  for  light  projectiles  inelastic  energy  losses  must  be  considered. 
The  stopping  powers  can  be  used  to  obtain  the  range  of  the  ions  as  follows: 


p(e)  = J 
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It  must  be  noted  that  this  equation  is  valid  only  when  the  initial  energy  of  the  ion 
is  much  greater  than  the  energy  transmitted  in  the  collision. 

2.2.1  Elastic  Collisions 

In  the  LSS  model  we  assume  a particular  interatomic  potential  between  the  ions 
and  atoms  in  the  solid.  Only  two  atoms  are  involved  in  a collision  according  to  the 
stopping  theory  and  this  is  called  the  Binary  Collision  Approximation  (BCA).  The  BCA 
is  justified  for  high  energies  when  the  two  colliding  electrons  are  at  a separation  much 
smaller  than  the  interatomic  distance  in  the  solid  and  the  probability  of  three  or  more 
particles  colliding  is  very  small.  Moreover,  the  collision  time  is  much  smaller  than  the 
vibrational  period  of  the  atoms.  The  interatomic  potential  will  depend  only  on  the 
distance  between  the  colliding  atoms.  With  reference  to  figure  2.2,  the  relationship 
between  the  transmitted  energy,  T and  the  scattering  angle  ®c 
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where  E is  the  initial  energy  of  the  projectile. 


2.2.2  Inelastic  Collisions 

For  the  case  of  low  velocities  and  amorphous  materials  it  can  be  assumed  that  the 
inelastic  energy  loss  is  proportional  to  the  ion  velocity.  From  the  work  of  Lindhard  and 
Scharff  [Lindhard  61],  the  electronic  energy  loss  can  be  approximated  as 
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The  value  of  k decreases  with  ion  mass,  which  means  that  the  inelastic  energy  losses  are 
more  important  for  lighter  ions. 

For  low  velocities,  the  scattering  resulting  from  a collision  is  between  the  ion  and 
the  electron  gas.  At  higher  ion  velocity,  the  ion  has  a probability  to  lose  all  its  electrons 
and  the  Lindhard  approximation  is  not  valid.  In  such  a case  the  stopping  power  will  be 
proportional  to  v'2  in  accordance  with  the  Bohr  model  [Bohr  48],  For  the  case  of  even 
higher  velocities,  inelastic  stopping  power  calculations  are  best  done  using  the  Bethe  and 
Bloch  [Townsend  76]  approximation. 


Since  the  time  of  Newton,  scientific  theories  have  been  developed  in  a 
reductionistic  mode,  i.e.  a complex  system  is  reduced  to  one  or  more  simple  subsystems 
which  are,  in  turn,  analyzed.  At  a higher  level,  simulations  serve  as  an  alternative  to 
reductionism  since  simulations  allow  us  to  study  the  behavior  of  classes  of  systems  or 
subsystems.  Therefore,  while  in  reductionism  the  emphasis  is  on  structural  analysis, 
simulations  stress  on  behavioral  classification.  In  fig.2.3,  it  is  shown  that  simulation  has 
two  distinct  roles  to  play  in  scientific  research.  At  a higher  level,  simulation  can  be  an 
alternative  to  reductionism  whereas  at  a lower  level  simulation  becomes  a tool  in 
reductionism. 


2.3  Computer  Simulations  and  Theory 


Analytical  theories  treat  the  issue  of  ion  bombardment  by  using  statistical 
quantities  such  as  cross  sections.  In  these  theories,  details  of  individual  ion  bombardment 
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are  not  clear.  On  the  other  hand,  in  computer  simulations,  single  ions  are  dealt  with  and 
their  interaction  with  the  target  solid  is  examined.  Experimental  tools  such  as  scanning 
tunneling  microscopy  and  atomic  force  microscopy  corroborate  some  results  of 
simulations  such  as  molecular  dynamics  simulations.  Moreover,  computer  simulations 
provide  a means  of  measuring  quantities  not  accessible  to  experiments. 

However,  the  limitation  of  computer  simulation  and  theory  starts  when  data  is 
needed  for  a wide  range  of  ion-target  combinations.  In  this  respect  range  theory  provides 
qualitative  values  for  all  ion-target  combinations.  Computer  simulations  can  be  either 
binary  collision  approximation  (BCA)  techniques  or  classical  molecular  dynamics.  BCA 
methods  are  computationally  very  efficient  and  can  be  employed  to  measure  quantities 
such  as  sputter  yields,  ranges  or  damage  profiles.  However,  the  BCA  models  make  use  of 
simple  potentials  unlike  classical  molecular  dynamics,  which  employs  more  accurate 
potentials  fitted  to  specific  properties  of  the  material  being  studied.  This  enables  one  to 
understand,  at  a microscopic  level,  the  processes  involved  in  the  bombardment  of  an  ion 
on  a particular  target.  The  trade  off  is  that  the  molecular  dynamics  simulations  are 
computationally  expensive.  It  follows  that  sputtering  yields  and  ranges  are  not  easily 
computed  with  this  technique. 

2.3.1  Binary  Collision  Approximation 

The  first  binary  collision  approximation  (BCA)  simulation  was  reported  by 
Yoshida  in  1961  [Yoshida  61],  The  BCA  models  can  be  classified  into  the  binary  crystal 
(BC)  models  and  the  Monte  Carlo  (MC)  models.  The  BC  model  employs  a crystalline 
target  with  all  the  atoms  having  well-defined  initial  positions.  MARLOWE,  developed  by 
Robinson  and  Torrens  [Robinson  74],  is  a popular  BC  model.  Modifications  and 
improvements  to  this  model  such  as  UT-MARLOWE  have  been  done  more  recently. 
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Damage  and  dopant  profiles,  sputtering  yields  and  a wide  range  of  target-projectile 
combinations  have  been  simulated  using  UT-MARLOWE  [Tasch  95],  Monte  Carlo 
models  use  random  methods  to  locate  target  atoms  as  well  as  random  values  to  obtain 
collision  parameters.  In  this  category,  TRIM,  a program  developed  by  Biersack  and 
Haggmark  [Biersack  80]  and  updated  yearly  by  Ziegler  is  very  popular  among 
experimentalists  and  in  the  modeling  world. 

In  order  to  use  time  scales  and  temperatures  typically  used  in  semiconductor 
manufacturing,  programs  for  calculating  damage  production  are  coupled  with  a Monte 
Carlo  simulation  of  defect  diffusion.  This  simulation  is  carried  out  by  first  generating  a 
cascade  by  MARLOWE.  The  locations  of  vacancies  and  interstitials,  obtained  from  the 
set  of  sites  where  a displaced  atom  was  in  its  final  and  initial  positions  respectively,  are 
passed  on  to  the  Monte  Carlo  diffusion  simulator.  Vacancy  -interstitial  recombination, 
annihilation  at  the  surface,  clustering  of  like  defects,  re-emission  from  the  clusters  and 
interstitial  trapping  at  native  carbon  traps  are  achieved  by  giving  the  vacancies  and 
interstitial  random  jumps.  The  accuracy  of  the  MC  diffusion  simulations  depends  on  the 
accuracy  of  the  parameters  used.  It  should  be  noted  that  Binary  collision  models  have 
succeeded  in  generating  subcascades  for  high  incident  particle  energies  [Beeler  76], 
2.3.2  Molecular  Dynamics  Simulations 

Molecular  Dynamics  (MD)  simulations  integrate  the  equations  of  motion  for  all 
the  atoms  in  a computational  cell.  In  classical  MD  simulations  interaction  between  the 
particles  is  determined  by  empirical  interatomic  potentials.  Particle  trajectories  are 
calculated  by  integrating  the  Newton’s  equations  of  motion.  The  gradient  of  a 
conservative  potential  yields  the  force  as  a function  of  the  particle  coordinates.  The 
relevant  equations  are  listed  below: 
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Fx  = — W(r) 
d2r  F 


(2.8) 

(2.9) 


dt2  mi 

These  equations  can  be  numerically  integrated  to  give  the  particle  trajectories  using  time 
as  the  integration  step.  There  is  a trade-off  between  computation  efficiency  and  errors  in 
choosing  the  time  step.  Usually,  the  time  step  is  of  the  order  of  0. 5-1.0  femtoseconds  (fs). 

2.3.2. 1 Newtonian  Dynamics 

For  N spherical  molecules,  Newton’s  second  law  represents  3N  second  order 
ordinary  differential  equations  of  motion.  MD  simulations  require  a numerical  method  for 
solving  the  differential  equations  of  motion.  The  classical  tools  for  solving  initial  value 
problems  are  finite  difference  methods.  These  methods  replace  differentials  such  as  dx 
and  dt  with  finite  differences  Ax  and  At.  Therefore  they  replace  differential  equations 
with  finite  difference  equations.  Also,  over  a small  but  finite  time  interval  At,  these 
methods  assume  that  the  rate  or  some  known  function  of  the  rate  is  constant 

2.3.2.2  Periodic  Boundary  Conditions 

Typically  MD  simulations  are  done  on  systems  containing  several  hundred  or  a 
few  thousand  atoms.  In  simulations  where  surface  effects  are  not  of  interest,  they  are 
removed  using  periodic  conditions.  Periodic  boundary  conditions  are  implemented  in  a 
simulation  of  N atoms  confined  to  a volume  V,  by  assuming  that  volume  V is  only  a 
small  fraction  of  the  bulk  material.  This  volume  V,  called  the  primary  cell,  is 
representative  of  the  bulk  material  with  the  assumption  that  the  bulk  is  composed  of  the 
primary  cell  surrounded  by  exact  replicas  of  itself  called  image  cells.  Therefore  the 
primary  cell  is  imagined  to  be  periodically  replicated  in  all  directions  to  constitute  a 
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macroscopic  sample  of  the  substance  of  interest.  This  periodicity  extends  to  the  positions 
and  momenta  of  the  images  in  image  cells. 

2.3 .2.3  Predictor-Corrector  Algorithms 

Predictor-Corrector  algorithms  consist  of  three  steps:  prediction,  evaluation  and 
correction.  Specifically,  from  the  current  position  x(t)  and  velocity  v(t)  the  steps  are  as 
follows: 

1 . Predict  the  position  x(t+At)  and  velocity  v(t+At)  at  the  end  of  the  next  step. 

2.  Evaluate  the  forces  at  t+At  using  the  predicted  position. 

3.  Correct  the  predicted  positions  and  their  derivatives  using  the  discrepancy  between  the 
predicted  acceleration  and  that  given  by  the  evaluated  force  (using  the  potential  energy 
function). 

In  this  work,  the  MD  simulation  code  used  is  called  MDCASK.  This  code  is 
adaptable  to  work  stations.  The  classical  equations  of  motion  are  integrated  numerically 
using  a fourth  order  predictor-corrector  scheme.  The  Langevin  equation  of  motion  is 
applied  to  the  atoms  in  the  link  cells  adjacent  to  the  cell  boundaries  in  order  to  control 
crystal  temperature. 

2. 3 .2. 4 Interatomic  Potential 

The  interatomic  potential  contains  the  information  about  the  interaction  between 
the  particles.  The  first  potentials  to  be  used  were  the  Born-Mayer  potential  [Gibson  60] 
and  the  Lennard-Jones  potential.  More  sophisticated  potentials  have  since  been  used. 
Examples  are  the  embedded  atom  poential  (EAM)  by  Daw  and  Baskes  [Daw  84]  which 
has  proven  to  be  very  successful  in  cascade  simulation  in  metals  and  alloys  [Diaz  de  la 
Rubia  92].  Some  of  the  most  important  potentials  for  covalent  materials  are  the  Stillinger- 
Weber  potential  [Stillinger  85],  the  Tersoff  potential  [Tersoff  89]  and  more  recently,  the 
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environment-dependent  interatomic  potential  (EDIP)  [Bazant  91].  A discussion  of  these 
potentials  for  silicon  can  be  found  in  later  chapters. 

Prior  to  studies  of  ion-implantation  in  semiconductors,  radiation  damage  studies 
in  metals  were  carried  out  beginning  with  Vineyard  and  co-workers  at  Brookhaven 
[Vineyard  72].  Related  studies  in  collision  cascade  phenomena  as  a result  of  radiation 
damage  have  been  done  since  with  the  advent  of  powerful  computers. 

To  summarize,  MD  appears  to  be  the  most  powerful  tool  in  the  study  of  ion- 
implantation  of  ions  with  energies  between  eV  and  a few  keV  for  the  following  reasons: 
MD  avoids  the  approximations  which  are  found  in  analytical  theories.  MD  simulations 
incorporate  the  target  structure.  MD  can  be  used  to  study  the  defect  distribution  in  the 
cascade.  Using  empirical  potentials,  only  reasonable  computing  times  are  required  for 
simulations  with  millions  of  atoms 

In  this  work  only  classical  MD  has  been  used  and  the  code,  MDCASK,  which  is  used  to 
carry  out  the  simulations  will  be  discussed  in  the  next  chapter. 
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Fig.  2.1  Nuclear  and  Electronic  stopping  powers  as  a function  of  reduced  energy  8 


Fig.  2.2  Two  atom  scattering 
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Fig.  2.3  Hierarchy  of  scientific  methods  of  investigation 


CHAPTER  3 

COMPUTER  CODE  FOR  MD  SIMULATIONS:  MDCASK 
MDCASK  is  a code  originally  based  on  the  MOLDY  code  developed  at  Harwell 
laboratory  in  the  United  Kingdom  by  Finnis  and  co-workers  [Finnis  89],  MDCASK  can 
be  used  to  study  ion  implantation  in  the  energy  range  between  eV  and  a few  keV  in 
silicon,  silicon  carbide  and  metals  targets. 

MD  simulations  generally  require  an  algorithm  for  integrating  the  equations  of 
motion,  specific  boundary  conditions  and  an  interatomic  potential.  Integration  of 
equations  of  motion  is  carried  out  in  MDCASK  using  a fourth-order  predictor  corrector 
algorithm  [Berendsen  85]  which  is  explained  in  detail  in  the  next  section. 

There  are  as  many  as  30  interatomic  potentials  today  for  silicon.  In  this  work  the 
Stillinger  Weber  potential  is  primarily  used  although  some  work  has  been  done  with  the 
Tersoff  and  the  Environment  Dependent  Interatomic  Potential  (EDIP)  which  will  be 
presented  in  subsequent  chapters.  MDCASK  uses  periodic  boundary  conditions  whereby 
the  MD  computational  cell  is  repeated  to  simulate  an  infinite  solid  [Born  12].  This  is 
accomplished  through  the  minimum  image  conversion  method  [Allen  86],  In  a real  solid, 
if  implantation  is  carried  out,  the  initial  ion  energy  will  spread  out  in  the  solid  by 
conduction  resulting  in  thermal  equilibrium  with  the  same  initial  temperature.  On  the 
other  hand,  in  simulations,  upon  the  application  of  boundary  conditions  the  energy 
remains  inside  the  solid  increasing  its  temperature.  The  MDCASK  code  has  features  to 
achieve  constant  temperature  as  in  a real  solid  which  will  be  explained  in  this  chapter. 
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3.1  The  Predictor  Corrector  Algorithm 

The  technique  to  integrate  the  equations  of  motion  used  in  MDCASK  is  called  the 
fourth-order  predictor  corrector  scheme.  Assuming  that  the  classical  trajectory  of  the 
atoms  in  the  system  is  continuous,  a Taylor  expansion  can  be  used  to  calculate  the 
positions  (r),  velocity  (v),  acceleration  (a),  etc.  at  a time  t+8t,  given  the  value  of  these 
quantities  at  time  t.  The  expressions  for  the  dynamic  variables  would  be  predicted  as 
follows,  if  the  Taylor’s  series  is  truncated  in  the  third  derivative: 


rp(t  + St)  = r(t)  + 8tv(t)  -l-  — 5t2a(t)  + — 8t3b(t) 

2 6 

vp(t  + St)  = v(t)  + 5ta(t)  + — 8t2b(t) 

2 

ap(t  + St)  = a(t)  + 8tb(t) 
bp(t  + St)  = b(t) 


(3.1) 

(3.2) 

(3.3) 

(3.4) 


The  predicted  values  for  the  positions  lead  to  the  forces  at  time  t+8t  using  the  inter- 
atomic potential.  The  correct  accelerations,  ac(t+8t),  give  us  the  error  in  the  prediction 
step  as: 


Aa(t  + 8t)  = ac(t  + 8t)  - ap(t  + 8t) 


(3.5) 


This  error  is  incorporated  in  the  predicted  values,  thus  yielding  the  corrected  new 

positions  and  velocities  as 

rc(t  + St)  = r^t  + 8t)  + c0Aa(t  + 8t) 

vc(t  + 8t)  = vp(t  + 8t)  + c^Aalt  + 8t) 

ac(t  + St)  = ap(t  + St)  + c2Aa(t  + 8t) 

bc(t  + 8t)  = bp(t  + St)  + c3Aa(t  + 8t) 


(3.6) 

(3.7) 

(3.8) 

(3.9) 
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The  values  of  the  coefficients  Co,  Ci,  C2  and  C3  have  been  optimized  by  Gear  and  are  the 
following  for  a fourth  order  predictor  corrector  algorithm: 

Co  = 1/6,  Ci  = 5/6,  C2  = 1 and  C3  = 1/3 

The  finite  time-step  of  the  simulation  introduces  errors  in  the  solutions  of  the  classical 
equations  of  motion.  Therefore,  the  optimum  time  step  is  chosen  such  that  At«l/u 
where  t)  is  the  vibrational  frequency  of  the  atoms  in  the  solid.  In  MD  this  time  step  is  a 
femtosecond. 


3.2  Cell  Structures  and  Linked  Lists 

In  the  beginning  of  the  simulation,  a list  is  constructed  of  all  the  neighbors  of  each 
molecule.  In  order  to  facilitate  this,  the  cubic  simulation  box  is  divided  into  smaller  cells, 
the  dimensions  of  each  being  slightly  larger  than  the  cut  off  distance  of  the  interatomic 
potential.  For  example,  as  shown  in  fig.  3.1,  the  neighbors  of  the  molecules  in  the  two 
dimensional  cell  13  are  found  in  cells  7,  8,  9,  12,  13,  14,  17,  18,  19.  The  linked  lists 
method  can  be  adopted  while  using  the  cell  structure  described  above.  The  first  step 
involves  sorting  all  the  molecules  into  their  appropriate  cells.  The  identification  of  one  of 
the  molecules  in  each  of  the  cells  is  stored  in  an  array  while  executing  the  program.  This 
head-of-chain  molecule  number  is  used  to  identify  the  element  of  a linked  list  array 
which  contains  the  number  of  the  next  molecule  in  that  cell.  This  linked  list  array  element 
for  this  molecule,  in  turn,  indicates  the  next  molecule  in  the  cell  and  this  process  is 
continued  until  an  element  of  the  linked  list  array  becomes  zero.  At  this  juncture,  it 
becomes  clear  that  there  are  no  more  molecules  in  that  cell  and  the  code  moves  on  to  the 
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head-of-chain  molecule  of  the  next  cell.  The  MD  program  calculates  the  forces  by 
looping  over  all  cells. 


Here,  the  Stillinger-Weber  interatomic  potential  is  discussed.  The  two  other 
interatomic  potentials  used  in  this  work  namely,  the  Tersoff  potential  and  the 
environment-dependent  interatomic  potential  (EDIP)  will  be  described  in  another  section. 
Any  interatomic  potential  must  have  terms  denoting  interactions  between  one  body,  two 
body  and  three  body  terms  such  as 


Ever  since  the  Stillinger-Weber  (SW)  potential  was  developed  by  Stillinger  and 
Weber  in  1985  [Stillinger  85],  it  has  been  used  extensively  in  numerous  computer 
simulations.  This  potential  describes  very  well  the  solid  and  liquid  states  in  silicon.  Since 
the  Si  crystal  consists  of  strong  and  directional  bonds,  it  would  be  appropriate  to  develop 
a potential  energy  function  comprising  pair  (u2)and  triplet  (d3)  terms  and  this  is  precisely 
the  form  adopted  by  Stillinger  and  Weber.  Upon  introducing  an  energy  unit  eand  a length 
unit  G,  the  pair  and  triplet  terms  can  be  written  as 


3.3  The  Stillinger-Weber  Interatomic  Potential 


^ = X W(i)  + X j)  + X i k> 


(3.10) 


(3.11) 


t)2(rij)  = £f2(ri;j  / G) 

v3(rif  rjf  rk)  = ef3(rL  / a,  r,  / G,  rk  / a) 


(3.12) 
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Here  8 is  chosen  to  give  f2  depth  -1,  and  a is  chosen  to  make  f2(2l/6)  vanish.  The  reduced 
pair  potential  used  by  the  authors  of  the  SW  potential  was  selected  from  the  following 
five-parameter  family: 


f2(p) 


{ 


A(Bp'p-p'q)exp[(p-a)‘  ] 


0 


r < a 
r>  a 


(3.13) 


where  p = r/a.  This  form  enables  us  to  cut  off  the  function  at  p = a,  with  the  derivative 
being  continuous.  The  three  body  term  selected  by  the  authors  with  the  same  cutoff 
advantage  is  as  folows: 


^(PbPjiPk)  h(pjj,pjk,0jik)  + h(pjj,  pjk,  0jjk)  + h(pkj,  pkj,  0jkj)  (3.14) 

In  this  equation,  0jlk  is  the  angle  between  pj  and  pk  subtended  at  vertex  i.  The  h functions 
take  the  form 


h(pij,pik,0jik)  = ^expfyCp.j-  a)'1  + [y(pik-  a)'1]  (cos0Jlk+l/3)2  (3.15) 

The  most  satisfactory  parameter  set  to  be  fitted  to  this  potential  is  as  follows: 

A = 7.049556277,  B = 0.602245584, 

P=4,  q = 0,  a = 1.8, 

A.  = 21.0,y=  1-2 


The  parameters  listed  above  were  obtained  such  that  the  melting  point  and  the 
liquid  structure  computed  by  MD  simulations  for  silicon  conform  to  experimental  values. 
This  potential  gives  the  melting  point  of  silicon  as  169  IK  within  an  error  of  20K.  which 
is  very  close  to  the  experimental  value  of  1683K.  In  addition  to  this,  the  diamond 
structure  is  proposed  to  be  the  most  stable  at  low  pressure.  The  cut  off  for  this  potential  is 
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found  to  be  at  a distance  3.77  A which  is  between  the  first  and  second  nearest  neighbor 
distance  in  diamond  silicon. 

Compared  to  experimental  values,  SW  yields  a value  for  the  bulk  modulus  very 
close  to  the  experimental  results.  As  will  be  evident  in  subsequent  chapters,  SW  yields 
good  structure  and  energetics  information  of  point  defects.  Also,  SW  matches  the  first 
principles  calculations  in  producing  the  (100)  Si  surface, but,  fails  to  model  of  the  the 
several  reconstructions  of  the  (111)  surface.  Since  the  SW  potential  consists  of  only  first 
and  second  nearest  neighbor  terms,  it  cannot  distinguish  between  the  diamond  phase  and 
the  hexagonal  diamond  phase.  By  far,  the  SW  potential  is  the  most  widely  used  potential 
in  simulations  of  silicon. 

3.3  The  Universal  Potential 

To  give  a better  treatment  of  interactions  between  atoms  at  short  distances  or 
having  higher  energies  the  repulsive  part  of  the  two-body  term  of  the  SW  potential  has 
been  modified  to  link  it  to  the  Universal  potential.  Although  there  is  no  sharp 
demarcation  for  the  validity  of  the  two  potentials,  this  range  is  chosen  to  be  between  a 
few  eV  for  the  potential.  In  the  MDCASK  code,  the  universal  potential  and  the  SW 

potential  are  connected  between  the  r values  in  the  range  n = 1 .2 1 07  A and  r2  = 1 .8456 

A.  These  r values  correspond  to  a potential  energy  = 24  eV  and  the  pair  part  for  the  SW 
equal  to  0 eV  for  n and  r2  respectively.  A spline  is  used  to  connect  the  two  potentials 
between  these  two  points. 

For  the  total  distance  range  the  value  of  the  pair  potentials  is  given  by 


V(r)  = 


(3.16) 


a0  + at  r + a2  r2+a3  r3ri  < r < 
e f2(r/a) 


T2 

r > i2 


where  a0  = 335.0594  eV,  a,  = -492.3905  eV  A a2  = 244.5280  eV  A ‘2,  and  a3  = - 
41.2347  eV  A ‘3. 

3.4  Constant-Temperature  MD  and  the  Canonical  Ensemble 

In  a MD  simulation,  the  positions  and  momenta  of  all  the  particles  form 
the  phase  space.  The  phase  space  points  constitiute  the  MD  ensemble.  An  average  over 
phase  space  gives  any  desired  property  of  the  system.  The  efficacy  of  the  MD  simulations 
in  producing  valid  results  depends  on  the  assumption,  from  a statistical  mechanics 
perspective  that  properties  or  observables  are  averages  over  some  region  of  sample  space. 
The  transition  from  a single  MD  system  to  an  ensemble  is  based  on  the  quasiergodic 
hypothesis  which  states:  ensemble  averages  are  equal  to  time  averages,  if  the  average  is 
done  for  a period  of  time  sufficient  to  make  the  average  independent  of  the  averaging 
time  [Beeler  76,  Heermann  86].  It  follows  that  in  order  to  build  statistical  ensembles  in 
MD  the  system  must  reach  an  equilibrium.  When  the  mean  and  variance  of  a dynamical 
variable  become  independent  of  time  equilibrium  is  attained. 

The  thermodynamic  state  of  a system  is  determined  from  three  parameters  such  as 
the  number  of  particles,  N,  the  temperature  T and  pressure,  P.  The  microcanonical 
ensemble  is  the  regular  ensemble  in  an  MD  simulation.  In  this  ensemble,  the  number  of 
particles,  N,  the  volume,  V and  the  total  energy  E are  constant.  In  this  ensemble,  the 
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problem  is  totally  deterministic  since  no  stochastic  elements  are  included  excepting  for 
the  elements  related  to  the  machine  round-off  errors.  This  is  an  isolated  system. 

In  many  cases,  however,  it  is  interesting  to  study  an  ensemble  with  N,  V and  T 
constant  or  N,P  and  H constant.  Constant  temperature  calculations  are  possible  in  the 
MDCASK  code  and  two  methods  of  achieving  this  are  highlighted  here. 

The  NVT  ensemble  is  known  as  the  canonical  ensemble.  In  this  ensemble,  the 
equations  of  motion  are  altered  to  couple  the  system  to  a thermal  bath.  In  this  case  only 
the  kinetic  energy  is  conserved.  The  equations  of  motion  are  constrained  by  a stochastic 
force  which  simulates  collisions  with  virtual  particles.  The  other  means  of  maintaining 
the  temperature  constant  is  by  fixing  the  kinetic  energy  of  the  system.  This  can  be 
achieved  by  rescaling  the  velocities  at  each  time  step  by  a factor  (T/t)1/2  where  x is  the 
current  kinetic  temperature  and  T is  the  desired  temperature  [Allen  87], 

Another  means  of  achieving  the  canonical  ensemble  is  by  using  Brownian 
dynamics.  In  this  technique,  the  system  particles  are  interacting  with  a viscous  medium. 
A stochastic  force  acting  on  every  particle  in  the  computational  cell  represents  the 
thermal  bath.  The  particles,  in  this  case,  are  said  to  obey  a Langevin  equation  of  motion 
in  a field  force  such  as 

m ^ = FA  + R±(t)  - K (3.17) 

dt 

where  (3vj  is  the  friction  force  acting  on  the  particle,  with  [3  being  the  viscous  constant 
and  Rj(t)  is  a random  force  that  simulates  the  heat  of  the  particle  at  the  necessary  bath 
temperature.  The  friction  force  removes  the  excess  of  energy  and  the  random  force 
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maintains  the  system  at  the  required  temperature.  The  random  force  must  conform  to 
certain  rules  such  as  the  following: 

It  has  to  be  zero  on  the  average  over  long  times,  i.e. 

< R(t)  > = 0 

The  random  force  must  be  uncorrelated  at  two  different  times  and  taking  into 
account  the  equipartition  theorem  must  obey 

<R(t)R(0)>=  2(3kT8(t)  (3.18) 

Moreover,  in  order  to  get  a Maxwellian  velocity  distribution  the  random  force 
must  be  Gaussian  distributed  as  follows: 

1 , , (3.19) 

P(R)  = t"  exp  (— R2  / 2 < R2  » 

V27t  < R2  > 

The  variance  <R2>  is  obtained  given  the  fact  that  the  correlation  time  for  the 
force  is  8(t)  as  follows: 

<R2>  = (3kBT/8t  (3.20) 

In  order  to  satisfy  these  conditions,  a random  number  must  be  generated  with 


mean  value  zero  and  variance  <R2>. 
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Fig.  3.1  The  cell  method  in  two  dimensions,  (a)  The  division  of  the  central  box  into  MxM 
cells  where  M=5.  (b)  Close-up  view  of  cells  1 and  2,  indicating  the  molecules  and  the 
link-list  structure 


CHAPTER  4 

DIFFUSIVITY  OF  POINT  DEFECTS  IN  SILICON:  A COMPARATIVE  STUDY 
USING  THREE  DIFFERENT  INTERATOMIC  POTENTIALS 

4. 1 Overview  of  Empirical  Interatomic  Potentials 
Normally,  empirical  potentials  are  constructed  by  first  guessing  a functional  form, 
based  on  physical  intuition,  and  then  by  choosing  the  parameters  to  fit  ab  initio  energy 
data  . This  process  is  particularly  difficult  in  covalent  materials  due  to  the  various 
complex  quantum-mechanical  effects  such  as  chemical  bond  formation  and  rupture, 
hybridization,  metallization,  charge  transfer  and  bond  bending.  There  is  a large  number 
of  potentials  for  silicon  in  the  literature.  Although  the  successes  and  failures  of  the 
various  potentials  cannot  be  attributed  to  the  functional  form,  it  can  be  concluded  that  a 
good  functional  form  is  more  important  than  intricate  fitting  strategies. 

As  described  in  the  previous  chapter,  the  Stillinger-Weber  (SW)  potential  has 
only  eight  parameters  fitted  to  a few  experimental  properties  [Stillinger  85],  This 
potential  has  a pair  interaction  term  which  has  a deep  well  at  the  first  neighbor  distance  to 
account  for  the  fact  that  there  is  a strong  restoring  force  opposed  to  the  stretching  of  the 
hybrid  sp3  covalent  bonds.  The  three-body  term  is  represented  as  a product  of  a radial 
and  an  angular  function.  The  angular  function  has  a minimum  of  zero  at  the  tetrahedral 
angle  suggesting  that  there  is  an  angular  preference  of  sp3  bonds.  On  the  other  hand,  the 
radial  function  decreases  with  distance  in  order  to  mitigate  the  effect  of  the  angular  term 
upon  bond  stretching.  It  is  also  seen  that  although  the  various  terms  lose  their  relevance 
for  distortions  of  the  diamond  lattice  leading  to  the  elimination  of  sp3  hybridization,  the 
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SW  potential  happens  to  be  the  most  commonly  used  potential,  mainly  due  to  its 
simplicity  and  sound  physical  content. 

The  other  popular  empirical  potential  is  the  Tersoff  potential,  with  three  versions 
named  T1  [Tersoff  86],  T2  [Tersoff  88]  and  T3  [Tersoff  88],  The  T1  version  has  only  six 
adjustable  parameters  fitted  to  a small  bulk  poly  type  database.  Elastic  properties  were 
further  improved  in  subsequent  versions  by  including  seven  more  parameters.  The 
fundamental  difference  between  the  Tersoff  and  SW  potentials  is  that,  in  the  Tersoff 
potential,  the  strengths  of  individual  bonds  is  influenced  by  the  presence  of  surrounding 
atoms.  The  energy  is  the  sum  of  a repulsive  pair  term,  0R(r)  and  an  attractive  interaction 
term  p(Q  <E>a(t)  which  depends  on  the  local  bonding  environment,  denoted  by  a scalar 
quantity 


Here,  the  function  p(Q  represents  the  Pauling  bond  order.  The  three-body  term 
has  a form  similar  to  the  one  found  in  the  SW  potential,  excepting  for  the  fact  that  the 
angular  function  may  not  have  a minimum  at  the  tetrahedral  angle.  It  is  due  to  the  nature 
of  the  angular  forces  that  the  Tersoff  potential  does  not  have  an  overall  advantage  over 
the  SW  potential.  Most  of  the  potentials  in  the  literature  have  formats  similar  to  the  SW 
or  tersoff  formats.  However,  a number  of  these  potentials  have  functional  forms  have 
restricted  validity  . As  an  example,  the  potential  developed  by  Pearson  et  al  [Pearson  84], 
has  no  physical  motivation  and  is  simply  an  example  in  fitting.  The  authors  of  this 
potential  use  Lennard-Jones  two-body  terms  and  Axilrod-Teller  three  body  terms  which 


E = + PtCi^Ri,)] 


(4.1) 

(4.2) 


ij 


Ci3  = X V3(Rij' 
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has  no  relevance  to  covalent  materials.  As  another  example,  four-body  terms  are  included 
in  some  other  potentials  [Mistriotis  89]  with  no  significant  advantage. 

The  newly  developed  environment-dependent  interatomic  potential  [EDIP] 
[Bazant  97]  was  developed  keeping  these  factors  in  mind.  The  premise  of  the  authors  was 
that  an  improvement  over  the  SW  and  Tersoff  potentials  could  be  achieved  by  replacing 
simple  functional  form  with  more  flexible  forms  and  also  by  carrying  out  more  involved 
fitting  schemes.  Another  consideration  was  that  it  may  be  impractical  to  attempt  a 
simultaneous  fit  of  all  important  atomic  structures  (bulk  crystalline,  amorphous  and  liquid 
phases,  surfaces  and  clusters)  because  different  features  of  bonding  are  behind  different 
types  of  structures.  The  authors  also  believed  that  the  local  bonding  environment  is 
important  while  considering  surface  and  bulk  bonding  preferences.  The  goal  of  the 
authors  was  to  obtain  the  best  possible  description  of  condensed  phases  and  defects  with 
a theoretically  justified  functional  form 

The  EDIP  potential  uses  a theoretically  motivated  functional  form 
which  stresses  on  chemical  and  physical  trends.  This  functional  form  is  fitted  to  a small 
ab  initio  database  with  only  thirteen  parameters.  The  authors  of  this  potential  claim  that 
this  is  a significant  improvement  over  existing  models  and  local  structures  and  extended 
defects.  The  authors  also  claim  that  this  potential  describes  point  defects  in  the  bulk,  the 
concerted  exchange  path  for  self-diffusion  and  elastic  properties  of  bulk  silicon.  The 
other  claims  are  that  the  potential  predicts  core  structures  of  partial  dislocations  in  the 
glide  set  { 1 1 1 } which  agrees  very  well  with  ab  initio  results. 

The  functional  form  of  the  EDIP  potential  is  as  follows: 
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The  configuration  energy  is  expressed  as  a sum  over  single-atom  energies. 


E = XiEi  ’ each  containing  two  and  three-body  terms, 

Ei  = X V2<Ri:'  Zi)  + X X V3<Ri:'  Rik'  Zi) 


(4.3) 


where  V2(Rijf  Zi ) is  representative  of  pairwise  bonds  between  atoms  i and  j and 

V3(Rijf  Rik'  zi ) represents  angular  forces  between  i,j  and  k centered  at  atom  i.  These  two 

interactions  depend  on  the  local  environment  of  atom  i through  its  effective  coordination 
number  given  as 


In  this  equation  f(Rim)  is  a cutoff  function  indicating  the  contribution  of  neighbor  m to  the 
coordination  of  atom  i in  terms  of  the  separation  Rim.This  function  is  unity  for  r<c, 
dropping  gradually  from  1 to  0 between  c and  a and  is  0 for  r>a.  This  is  represented  as: 


where  x = (r-c)/(a-c).  For  r<c  a neighbor  of  atom  i is  considered  a full  neighbor,  whereas 
neighbors  between  c and  a only  partially  contribute  to  Z,  The  cutoffs  are  chosen  to 
reproduce  the  important  crystal  structure  coordinations  such  as 
Zj  = 4 for  the  diamond  lattice. 

The  two-body  term  includes  interactions  which  are  attractive  and  repulsive: 


zi  = Xf(RiJ 


(4.4) 


exp(a/(l-x'3)) 

0 


if  r<c, 
if  c<r<a 
if  r>a 


(4.5) 


V2(r,Z)  = A[(B/r)p-p(Z)]  exp(o/(r-a) 


(4.6) 


This  function  has  continuous  derivatives  and  goes  to  0 at  r = a It  should  be  noted  that 
although  this  choice  is  very  similar  to  the  SW  two-body  term,  the  unique  property  of  this 
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function  is  that  the  bond  strength  changes  according  to  the  local  atomic  environment. 
Theoretical  calculations  show  a weakening  of  the  attractive  interaction  and  therefore  an 
increase  in  bond  length  for  increasing  coordination  [Bazant  97],  A Gaussian  function 

demonstrates  this  theoretical  dependence  such  as 

n,2  (4.7) 

P(Z)  = e-pz 

The  three-body  term  consists  of  radial  and  angular  terms  and  is  as  follows: 


where  lyk  - cosByk  and  the  radial  function  is  selected  to  have  the  SW  form, 
g(r)  = exp(y/(r-a))  (49) 

The  angular  function  which  depends  strongly  on  the  local  coordination  by  means  oi  iwo 
functions  x(Z)  and  w(Z)  that  determine  the  equilibrium  angle  and  the  strength  of 
interaction  respectively.  By  theoretical  considerations  the  following  form  is  proposed: 


The  three  body  angular  function  becomes  flatter  as  the  coordination  increases. 
This  is  consistent  with  the  desire  to  have  a weakening  of  the  angular  forces  with 
increasing  coordination.  Therefore,  to  summarize,  this  version  of  EDIP  has  13  adjustable 
parameters  for  silicon  which  are:  A,  B,  p,  p,  c,  a,  c,  X,  p,  y,  p,  a and  Q 0.  One  of  the  other 
features  addressed  by  the  authors  is  that  of  computational  efficiency.  The  environment 
terms  in  the  two-  and  three-  body  terms  necessitate  extra  iteration  loops  in  the  force 
calculation.  Although  the  three-body  interaction  requires  a four-body  loop  leading  to  a 
more  expensive  force  evaluation  compared  to  the  SW  potential,  this  four-body  loop  needs 
to  be  performed  only  when  the  atoms  are  in  the  range  c<r<a.  Therefore,  the  authors 


(4.8) 


h(l,  Z)  = X 1 ~ e1 


) + T)Q(Z)  (1  + x(Z)  )2] 


with  w(Z)-2  = Q(Z)  = Q0  e'MZ 


(4.11) 
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claim,  that  increased  sophistication  and  realism  can  be  realized  with  minor  computational 
expense  in  the  EDIP  potential  compared  to  the  SW  potential. 


4.2  Diffusivities  of  Point  Defects  and  Self-Diffusion  in  Silicon 
It  has  been  observed  experimentally  that  the  silicon  self-diffusion  coefficient  shows 
Arrhenius  behavior  over  a wide  range  of  temperatures, 

D(T)  = D0exp(-Ea  / kBT)  ^ { 

where  the  activation  energy  Ea  is  in  the  range  4-5  eV.  The  preexponential  factor  D0  is 
much  greater  than  the  corresponding  values  in  typical  metals.  Elowever,  the  relative 
contribution  from  different  mechanisms,  such  as,  from  vacancies  and  interstitials  has  not 
been  resolved  yet.  Gosele  [Gosele  96]  has  recently  reported  that  the  experimental  data  for 
interstitial  (I)  and  vacancy  (V)  diffusion  are  as  follows: 

Cjdj  = 914  exp  (-4 . 84  / kBT)  cm2/sec  (4.13) 

and 

C*vdv  = 0 . 6 exp  (-4 . 03  / kBT)  cm2/sec  (4. 1 4) 

where  d,  and  dv  are  single  defect  diffusivities  and  C*  is  the  equilibrium  defect 
concentration  normalized  to  5 X 1022  cm'3.  It  has  also  been  found  [Gosele  96]  that  C]d;  is 
a more  reliable  and  accurate  quantity-  compared  to  the  vacancy  quantity  -since  consistent 
values  are  found  using  two  different  experimental  measurements.  Moreover,  although 
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Cjdj  is  well  established,  the  individual  values  of  the  terms  d and  C*  are  unknown  both 
for  interstitials  and  vacancies. 

In  this  work  the  diffusivities  di,v  are  calculated  theoretically  using  Molecular 
Dynamics  simulations  using  three  different  interatomic  potentials,  namely,  the  SW, 
Tersoff  and  EDIP  potentials.  The  values  of  the  migration  energy  barriers  of  the  interstitial 
and  the  vacancy  obtained  using  the  three  different  potentials  determine  the  validity  of  the 
potentials  for  molecular  dynamics  simulations. 

MD  calculates  the  atomic  trajectories  of  a group  of  interacting  atoms.  Since  the 
time  step  for  simulation  is  rather  low  (10'15  sec),  the  time  duration  of  simulation  is  rather 
low.  In  this  study  the  diffusion  is  simulated  for  1X1  O'9  sec.  The  diffusion  coefficients  are 
calculated  by  following  the  displacements  of  all  the  atoms  in  the  system.  The  accuracy  of 
the  results  is  increased  by  using  high  temperature  simulations.  The  computation  for  single 
defects  and  small  clusters  was  done  using  a computational  cell  consisting  of  216  atoms. 

A larger  computational  cell  of  1000  atoms  used  by  M.  J.  Caturla  at  Lawrence  Livermore 
National  Laboratory  yielded  results  almost  identical  to  the  results  from  the  216  atom  cell 
used  in  this  work. 

As  mentioned  earlier,  the  diffusion  coefficients  of  single  point  defects  and  small 
clusters  were  obtained  not  by  following  the  migration  of  the  defect  itself  but  from  the 
displacements  of  all  the  atoms  in  the  system.  As  an  example,  the  diffusion  of  a single 
vacancy  is  made  possible  by  the  nearest-neighbor  displacements  of  lattice  atoms  as  they 
are  transported  to  neighboring  vacant  sites.  In  the  absence  of  defects,  none  of  the  atoms  is 
moved  to  a neighboring  site.  The  diffusion  coefficients  of  the  point  defects  can  hence  be 
calculated  from  mean  squared  displacements  over  all  atoms, 
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D = ^[r^t)  - rt (0)  ] 2 / 6t  (4 A5) 

i 

This  equation  is  valid  for  larger  elapsed  time  t.  However,  due  to  the  large  number  of 
atoms  in  the  computational  cell  and  the  short  simulation  times,  the  vibrational 
displacements  of  atoms  around  their  lattice  sites  introduces  errors  in  the  computation. 
Specifically,  there  could  be  an  initial  steep  rise  in  the  plot  of  mean  squared  displacement 
as  a function  of  time  due  to  the  vibrational  displacements  of  atoms  around  their  lattice 
sites.  This  possible  error  can  be  eliminated  by  plotting  the  mean  squared  displacement, 
Xt  - ri  (0)  ] 2 as  a function  of  time,  and  determining  the  slope  of  the  curve  after  the 

i 

initial  steep  slope  of  the  curve,  due  to  the  vibrational  motion,  has  saturated  [Gilmer  95], 

A precaution  should  be  taken  about  the  orbital  motion  of  large  defect  clusters 
around  the  cluster’s  center  of  mass  while  calculating  the  diffusivity  of  the  clusters. 
Although  for  long  simulation  times  this  should  not  create  any  problem  since  the 
displacements  due  to  orbital  motion  saturates,  for  the  short  time  scales  of  MD  simulations 
the  orbital  motion  does  contribute  to  the  diffusivity  values  as  they  take  a longer  time  to 
saturate  than  the  vibrational  motion  of  atoms.  Fig.  4 1 illustrates  the  technique  of 
determining  the  diffusivity  from  the  mean  squared  displacement  versus  time  curve.  As  is 
evident  from  this  plot  there  is  a linear  relationship  between  the  mean  squared 
displacement  and  the  time  variable.  This  plot  is  for  a fixed  temperature.  These 
calculations  are  repeated  for  various  temperatures.  From  a semi-log  plot  of  the  diffusivity 
versus  1/kT,  where  k is  the  Boltzman’s  constant  and  T is  the  temperature,  the  migration 
energy  barrier  is  extracted.  This  method  is  adopted  for  calculating  the  migration  energy 
barriers  using  the  SW,  Tersoff  and  EDIP  potentials.  Fig.  4.2  is  a plot  of  temperature 
dependent  diffusivity  of  the  vacancy,  the  interstitial  and  the  di-interstitial  calculated  using 
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the  SW  potential.  Fig.  4.3  and  fig.  4.4  are  similar  plots  for  calculations  done  using  the 
Tersoff  and  EDIP  potentials  respectively.  Fig.  4.5  contains  Arrhenius  plots  of 
diffusivities  of  vacancy  for  the  three  different  potentials  in  order  to  show  a comparison 
between  the  potentials.  Figs.  4.6  and  4.7  are  similar  comparative  plots  of  interstitial  and 
di-interstitial  diffusivities  respectively  for  the  three  different  potentials. 

Table  4.1  contains  the  migration  energy  barriers  for  the  vacancy,  interstitial  and 
the  di-interstitial  for  calculation  done  using  the  S W,  Tersoff  and  the  EDIP  potentials. 
From  this  data  it  is  evident  that  the  EDIP  potential  gives  an  interstitial  migration  energy 
value  which  is  vastly  lower  (0.33  eV))  than  that  obtained  from  the  SW  (0.9  eV)  and 
tersoff  (1.1  eV)  potentials. 

First-principles  methods  also  yield  a value  for  the  interstitial  migration  energy 
close  to  that  obtained  from  the  SW  and  tersoff  potentials  [Blochl  93],  Blochl  et  al. 
obtained  the  self  diffusion  coefficient  by  combining  concentrations  and  diffusivities  . The 
concentrations  were  obtained  by  determining  the  formation  energies  using  ab  initio 
techniques  and  formation  entropies  from  the  local  harmonic  approximation  [Najafabadi 
89],  In  this  approximation,  for  a given  perfect  crystal  structure  at  volume  V and 
temperature  T,  one  can  determine  the  three  vibrational  frequencies  of  an  atom  by 
diagonalizing  the  local  dynamical  matrix  for  each  independent  atom  in  the  unit  cell. 

These  local  frequencies  yield  the  Helmoltz  free  energy  from  which  all  other 
thermodynamical  quantities  are  determined.  They  also  found  that  for  a migration  energy 
around  1 . 1 eV  for  the  self  interstitial  the  self  diffusion  values  are  in  the  range  of 
experimental  results.  This  therefore  suggests  that  the  migration  energy  values  obtained 
using  SW  and  tersoff  potentials  conform  to  first  principles  estimates. 
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The  authors  of  the  EDIP  potential  agreed  that  there  is  a flaw  in  this  potential. 
They  calculated  the  fully  relaxed  activation  pathway  between  the  tetrahedral  and 
hexagonal  sites  with  EDIP.  From  these  calculations  they  noticed  that  the  interstitial  path 
is  unrealistic  in  many  ways.  Mainly,  the  interstitial  path  has  an  artificial  global  minimum 


Table  4.1  Migration  energy  barriers,  of  vacancy  and  interstitial  calculated  for  the  SW, 
Tersoff  and  EDIP  potentials,  from  this  work. 


Potential 

Em  (Vacancy)  (eV) 

EM(Interstitial)  (eV) 

SW 

0.436 

0.926 

Tersoff 

0.383 

1.1 

EDIP 

0.203 

0.334 

between  the  tetrahedral  and  hexagonal  sites  and  perhaps  also  a low  activation  barrier.  The 
authors  attributed  this  flaw  to  a careless  error  that  had  been  introduced  into  the  fitting 
database.  A certain  intermediate  configuration  between  the  tetrahedral  and  hexagonal  was 
fitted  to  the  LDA  energy  of  the  hexagonal  minimum.  This  in  turn  decreased  the  true 
minimum  of  energy  and  also  apparently  created  an  invalid  global  minimum.  This  lack  of 
validity  of  the  EDIP  potential  means  that  this  potential  cannot  be  successfully  used  in 
studying  clustering  phenomena. 

In  order  to  further  investigate  the  validity  of  the  EDIP  potential  it  would  be 
necessary  to  calculate  the  self-diffusion  of  the  interstitial  and  the  vacancy.  A comparison 
of  the  temperature,  at  which  there  is  a crossover  between  the  interstitial-dominated  to  the 
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vacancy-dominated  self-diffusion  regime,  for  theoretical  calculations  using  the  EDIP 
potential  and  experimental  observations  would  further  indicate  the  validity  of  the  EDIP 
potential. 

A theoretical  estimate  of  the  self-diffusion  coefficients  requires  three  different 
calculations:  (a)  diffusivity  prefactor  d0  and  migration  energy  EM  (b)  formation  energy  Ef 
and  (c)  formation  entropy  Sf.  d0  and  EM  were  obtained  using  MD  with  the  EDIP  potential. 
The  self-diffusion  calculations  will  be  complete  only  when  equilibrium  concentrations 
are  computed.  The  formation  entropy  (Sf)  and  energy  (Ef)  of  the  individual  defects  are 
related  to  the  normalized  equilibrium  concentration  of  the  interstitial  and  vacancy 
through  the  equation 

C*i,v  = exp[-  (Ef-TSf)/kBT)  (4.16) 

The  calculation  of  formation  energy  is  described  in  the  following  section. 

4,3  Formation  Energy  Calculations  of  Point  Defects  and  the  Di-Interstitial 

The  formation  energy  of  the  point  defects  and  the  di-interstitial  was  obtained 
using  the  SW  and  the  EDIP  potentials.  For  self-diffusion  calculations  the  EDIP  potential 
was  used.  For  formation  energy  calculations  a computational  cell  of  1000  atoms  was 
used.  The  stable  structures  of  the  point  defects  were  obtained  first  by  annealing  the 
system  containing  the  defects  at  arbitrary  positions  at  1000K  for  50  picoseconds.  This 
tempeature  would  be  sufficient  to  produce  several  diffusion  hops  of  the  defects. 
Subsequently,  the  system  is  cooled  to  0 K over  a period  of  50  picoseconds  gradually. 

This  technique  is  found  to  give  a unique  low-energy  defect  structure  independent  of  the 
initial  structure.  The  validity  of  this  approach  is  further  corroborated  from  the  formation 
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energies  obtained  which  were  found  to  be  having  the  same  values  as  those  reported  in  the 
literature  for  the  SW  potential  [Gilmer  95]. 

4.4  Self-Diffusion  Data  for  the  EDIP  Potential 
Table  4.2  shows  the  formation  energy  values  obtained  in  this  work  using  the  SW 
and  EDIP  potentials  as  well  as  the  values  obtained  from  the  literature  for  the  Tight 


Table  4.2  Formation  energies  of  point  defects  and  the  di-interstitial  obtained  by 
simulations 


Defect  type 

SW 

TBMD 

ab  initio 

EDIP 

Vac. 

2.64  eV 

3.97  eV 

3.65  eV 

3.23  eV 

Int. 

3.65  eV 

3.8  eV 

3.7  eV 

3.43  eV 

di-int. 

5.7  eV 

4.91  eV 

" 

4.97  eV 

The  data  labeled  SW  and  EDIP  are  obtained  in  this  work.  TBMD  values  were  obtained 
from  Tang  et  al.  [Tang  97],  ab  initio  data  are  from  Zhu  et  al.  [Zhu  96], 

Binding  scheme  and  ab  initio  schemes.  The  self-diffusion  data  is  then  computed  using  the 
migration  energy  barrier  values  and  the  formation  energy  values  for  the  case  of  the  EDIP 
potential.  In  the  present  work,  the  formation  entropy  is  not  calculated  using  MD 
simulations.  For  vacancies,  Blochl  et  al.  have  calculated  the  formation  entropy  by  means 
of  the  local  harmonic  approximation  and  have  obtained  a value  of  Sfv  = 9 kB  [Blochl  93], 
For  interstitials,  Tang  et  al.  have  fitted  C*|d,  to  the  experimental  data  point  at  1250  °C  and 
have  extracted  the  formation  entropy  to  be  1 1.2  kB  [Tang  97],  C*id,  is  given  as 
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C*idi  = d0  exp(Sfl/kB)  exp[-(EffEm)/kBT]  ^ j ^ 

Similarly  the  self-diffusion  of  vacancies  is  given  by 

C\dv  = d0  exp(Sfv/kB)  exp[-(Ef+Em)/kBT]  <4- 1 8) 

Fig  4.8  depicts  the  Arrhenius  plot  of  the  self-diffusion  data.  From  this  plot  it  is 
estimated  that  the  crossover  temperature  from  the  interstitial-dominated  to  vacancy- 
dominated  self-diffusion  regime  falls  at  about  717°C.  This  temperature  value  is  vastly 
different  from  the  experimental  finding  of  1050°C.  Therefore  this  is  an  additional  proof 
that  the  EDIP  potential  is  lacking  in  validity.  The  SW  potential  and  the  Tersoff  potential 
give  crossover  temperatures  of  1 152  °C  and  1082  °C  respectively  which  are  closer  to  the 
experimental  value  than  the  value  obtained  using  the  EDIP  potential. 

4.5  Summary 

In  this  chapter  an  overview  of  SW,  Tersoff  and  EDIP  interatomic  potentials  was 
presented.  The  results  of  MD  calculations  of  the  migration  energy  of  point  defects  and 
the  di-interstitial  were  presented.  These  results  show  that  the  EDIP  potential  yields  a 
migration  energy  barrier  for  the  interstitial  which  is  much  lower  than  that  obtained  from 
SW,  and  Tersoff.  The  SW  and  Tersoff  values  are  very  close  to  that  obtained  from  LDA 
calculations,  thereby  implying  that  there  is  a flaw  in  the  EDIP  potential.  The  formation 
energies  of  the  point  defects  were  also  obtained  by  MD  simulations  using  the  EDIP 
potential.  From  the  formation  energies  and  migration  energies  the  self-diffusion 
coefficient  was  calculated  for  the  interstitial  and  the  vacancy  over  a range  of 
temperatures.  It  was  found  that  the  crossover  temperature  for  the  interstitials  and  vacancy 
self-diffusion  coefficients  was  717  °C  whereas  the  experimental  value  is  1050  °C  proving 
again  that  the  EDIP  potential  is  flawed. 
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y = 1.5421  + 4.1385e+10x  R=  0.9958 


Mean  squared  displacement  versus  time 


Fig  4.1  Mean  squared  displacement  (msd)  as  a function  of  time  for  interstitial  migration 
using  the  SW  potential.  This  is  for  the  case  of  the  interstitial  at  1450  K 
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— * — Vacancy(SW) 

— ▼-  - Interstitial(SW) 

- ■ - di-interstitial/2 


dv  = 0.00125  * eA(-0.43/(kBT))  R=  0.99 

d,=  0.0297  * eA(-0.92/(kBT))  R=  0.98 

d|2  = 0.0001 1 * eA(-0.15/(kBT))  R=  0.96 
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Fig.  4.2  Arrhenius  plot  of  diffusivity  as  a function  of  l/kBT  using  the  SW  potential 
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— ▼-  - Interstitial 
— ± — vacancy 
— ■-  - di-interstitial/2 


d,  = 0.032  * eA(-1. lAkgT))  R=  0.98 

— dv  = 0.00028*  eA(-0.38/(kBT))  R=  0.99 
■ d|2=  0.01*  eA(-1.1/(kBT))  R=  0.98 


Fig.  4.3  Arrhenius  plot  of  diffusivity  as  a function  of  l/kBT  using  the  Tersoff  potential 
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Fig.  4.4  Arrhenius  plot  of  diffusivity  as  a function  of  l/kBT  using  the  EDIP  potential 
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▼ s w 

— ■ -Tersoff 
<►-  EDIP 


dv=  0.00125  * eA(-0.43/(kBT))  R=  0.99 

- dv=  0.00028  * eA(-0.38/(kBT))  R=  0.99 

dv  = 0.00013  * eA{-0.2/(kBT))  R=  0.96 


Fig.  4.5  Arrhenius  plots  showing  vacancy  diffusivity  for  SW,  Tersoff  and  EDIP 
potentials 
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-A  s W 
— ■-  - Tersoff 
- ▼-  - EDIP 


d ,=  0.0297  * eA(-0.92/(kgT))  R=  0.98 

d,  = 0.032*  eA(-1.1/(keT))  R=  0.98 


d , = 0.00089*  eA(-0.33/(kBT))  R=  0.99 


Fig.  4.6  Arrhenius  plots  showing  interstitial  diffusivity  for  SW,  Tersoff  and  EDIP 
potentials 


Diffusivity  (cm2/sec) 


58 


d)2=  0.00011  * eA(-0.15/(kBT))  R=  0.96 

d|2  = 0.01*  eA(-1.1/(kBT))  R=  0.98 

d|2  = 0.0052*  eA(-0.62/(kBT))  R=  0.99 
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Fig.  4.7  Arrhenius  plots  showing  di-interstitial  diffusivity  for  SW,  Tersoff  and  EDIP 
potentials 
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— ▼ — vacancy 
— ■ - interstitial 


D = 1 .09  * eA(-3.43/(kBT)) 

- D(  = 53.3  * T)) 


Fig.  4.8  Arrhenius  plot  of  self-diffusion  for  vacancy  and  interstitial  using  the  EDIP 
potential 


CHAPTER  5 

ENERGETICS  OF  SMALL  SELF-INTERSTITIAL  CLUSTERS 

5.1  Introduction 

It  is  now  well  known  that  self-interstitial  atoms  in  silicon,  generated  by  ion- 
implantation  play  a significant  role  in  self-diffusion  and  impurity  diffusion.  Due  to  their 
high  formation  energies  (3-5  eV)  [Fahey  89],  self-interstitials,  also  aggregate  to  form 
clusters  and  extended  defects  such  as  the  {311}  defects.  Self-interstitial  agglomerates, 
such  as  {3 1 1 } defects  have  been  seen  to  be  involved  in  the  phenomenon  of  dopant 
Transient  Enhanced  Diffusion  (TED).  This  occurs  by  the  injection  of  interstitials  through 
the  dissolution  of  {31 1 } defects,  which  were  initially  generated,  upon  annealing  [ Stolk 
97],  Moreover,  it  has  been  observed  that  the  absolute  enhancement  and  duration  of  boron 
TED,  which  is  instigated  by  interstitial  supersaturation,  are  dependent  also  on  the  thermal 
stability  of  small  interstitial  clusters  [Zhang  95], 

Although  there  have  been  numerous  studies  on  obtaining  structural  information  of 
{311}  defects  using  Transmission  Electron  Microscopy  (TEM)  analyses  [Eaglesham  94], 
as  well  as  theoretical  studies  [Kohyama  92],  there  is  very  little  information  about  the 
structure  and  energetics  of  nanometer  -size  interstitial  clusters.  In  this  chapter,  the 
formation  energy  results  of  small  clusters  (size  < 1 0 interstitials),  calculated  using  MD 
simulations  with  the  SW  interatomic  potential,  will  be  presented.  In  addition  to  this,  the 
binding  energy,  obtained  from  the  formation  energy  calculations  as  well  as  the  atomic 
models  of  some  small  clusters  will  be  presented. 
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5.2  Computational  Methodology 

Classical  MD  simulations,  using  the  MDCASK  code  with  the  SW  potential,  are 
used  in  this  work.  Because  of  the  very  small  time  steps  involved  in  these  simulations,  the 
total  simulation  time  is  limited.  In  this  work,  a supercell  of  1 000  atoms  is  used.  First,  the 
case  of  a perfect  lattice  is  run  by  heating  the  perfect  lattice  to  1 000K  for  5 nanoseconds 
and  then  cooling  the  system  gradually  to  0 K.  Next  the  required  number  of  interstitials  is 
introduced  at  arbitrary  locations.  As  before,  the  system  is  annealed  for  5 nanoseconds  at 
1000  K,  and  subsequently,  it  is  cooled  slowly  over  a time  span  of  5 nanoseconds  to  0 K. 
This  procedure,  of  simulated  annealing,  has  the  advantage  that  the  atoms  will  undergo 
several  diffusion  hops  before  attaining  the  most  stable  configuration  after  relaxation. 

This  ensures  that  if  there  is  a metastable  configuration  for  the  clusters,  we  can  get 
out  of  it  by  quenching  the  temperature  slowly.  The  method  used  to  do  this  quench  is 
called  Langevin  dynamics  [Allen  87],  This  scheme  involves  coupling  the  atoms  to  a 
thermal  reservoir,  so  that  the  temperature  remains  constant.  The  program  does  this  in  the 
following  way:  after  the  calculation  of  the  force  on  each  atom,  there  is  a term  added  to 
that  force,  which  is  like  a friction  force  and  sets  the  temperature  essentially  to  0 K.  Then, 
another  term  is  added,  which  contains  a random  number,  so  that  the  desired  temperature 
is  added  to  the  system.  This  procedure  is  done  for  every  time  step.  The  advantage  with 
this  method  is  that  we  can  choose  how  fast  the  quench  is  occurring.  In  the  code,  this  is 
done  in  a subroutine  called  tfix.f 

The  input  file  for  the  MDCASK  code  is  called  mold.in.  In  mold. in,  there  are 
several  variables  that  are  related  to  the  procedure  described  above. 

Langevin  dynamics: 


TEMPCONT  NLANGEVIN  TEMPFACT 
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ZTC  XYTC  NTCLAYERS 
BETAL 

TEMPCONT:  If  this  variable  is  assigned  a logical  value  True,  then  this  temperature 
control  is  going  to  be  used. 

NLANGEVIN:  Determines  the  number  of  steps  after  which  the  temperature  should  be 
rescaled. 

TEMPFACT:  This  variable  decides  by  how  much  the  temperature  is  going  to  be  rescaled. 
For  example,  0.8  is  80%. 

This  temperature  rescaling  can  be  carried  out  for  either  all  the  atoms  in  the  simulation  or 
for  only  some  of  them.  The  following  parameters  determine  this  choice: 

ZTC:  If  true,  then  the  atoms  in  the  Z direction  will  be  scaled. 

YXTC:  If  true,  then  atoms  in  the  x and  y direction  will  be  scaled. 

NTCLAYERS:  This  sets  the  number  of  layers  of  atoms  to  be  scaled.  If  the  box  is 
10x10x10  and  20  layers  are  selected,  it  means  that  all  the  atoms  in  the  simulation  will  be 
coupled  to  the  reservoir. 

Supposing  that  an  interstitial  is  introduced  into  the  system  and  the  temperature  of 
anneal  is  set  at  1000  K,  the  code  starts  with  TEMPRQ  = 1000  K.  If  NLANGEVIN  = 

1000,  every  1000  steps  the  temperature  is  lowered  by  0.8  * TEMP,  if  TEMPFACT  = 0.8 
and  if  the  total  number  of  steps  is  stipulated  to  be  50000,  then  the  code  runs  for  50000 
steps. 

From  these  simulation  results  the  formation  energy  of  the  interstitial  is  obtained. 
The  formation  energy  of  a cluster  is  defined  as 
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E(Formation)  = Potential  energy  of  the  system  with  the  clusters  - (N+n)/N  * Potential 
energy  of  the  perfect  system. 

where  N is  the  total  number  of  atoms  in  the  perfect  box,  and  n is  the  total  number  of 
interstitials.  These  potential  energies  are  obtained  from  two  output  fdes,  namely, 
mold. out  and  mdyn.ene 

In  mdyn.ene  we  get  the  difference  between  the  total  energy  of  the  system  after 
relaxation  and  the  energy  at  the  beginning.  The  potential  energy  of  the  system  with  the 
cluster  is  = The  total  energy  from  mold.out+the  energy  from  mdyn.ene.  There  are  two 
ways  to  obtain  the  formation  energy  of  clusters.  In  the  first  method,  an  original  form  of 
the  cluster  is  set  up  and  the  temperature  is  increased  and  the  system  is  relaxed.  However, 
in  this  technique,  different  configurations  will  have  to  be  tried  to  find  the  lowest  energy 
configuration  from  all  of  them. 

In  the  second  method,  the  interstitials  are  introduced  at  random  positions  and  the 
system  is  taken  to  a high  temperature  and  held  for  a sufficiently  long  time  for  the 
interstitials  to  come  together  and  for  the  cluster.  The  disadvantage  with  this  method  is 
that  the  simulation  time  is  very  high.  In  this  work  this  second  technique  is  adopted. 

5.3  Results  and  Discussion 

Figure  5.1  shows  the  plot  of  formation  energy  per  interstitial  versus  defect  size. 
Figure  5.2  shows  a plot  of  binding  energy  versus  defects  size.  From  the  formation  energy 
per  interstitial  plot  it  is  evident  that  there  is  a significant  drop  in  formation  energy  as  the 
cluster  size  increases  especially  for  the  first  few  sizes.  This  implies  that  there  is  a driving 
force  for  the  aggregation  of  individual  interstitials  to  form  clusters.  Unfortunately, 
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because  of  computational  cost,  the  formation  energy  data  could  not  be  obtained  for  a 
defect  size  greater  than  9 interstitials. 

The  binding  energy  plot  shows  an  oscillatory  trend  with  increase  in  cluster  size. 
These  binding  energies  are  used  as  inputs  in  Monte  Carlo  calculations.  The  binding 
energies  calculated  in  this  work  agree  well  with  binding  energy  values  in  the  literature  for 
cluster  size  up  to  5 interstitials  [Jaraiz  96],  Jaraiz  et  al.  have  tabulated  binding  energy 
values  for  cluster  size  up  to  5 intersititials  and  the  oscillations  are  seen  in  their  work  as 
well.  Their  comment  was  that  although  oscillations  are  seen  at  very  small  cluster  sizes  in 
the  binding  energy,  there  should  be  a convergence  of  the  binding  energy  to  a single  value 
(2.1  eV)  at  very  large  cluster  sizes. 

Figure  5.3  shows  the  atomic  models  of  a two-interstitial  cluster.  As  is  seen  from 
this  figure,  there  is  a strain  field  around  the  cluster  and  lattice  atoms,  in  the  vicinity  of  the 
cluster  are  displaced  from  their  original  locations.  Figure  5.4  and  fig.  5.5  show  the  atomic 
models  of  a three-interstitial  cluster  and  a seven-interstitial  cluster  respectively.  In  these 
models  too  the  strain  field  around  the  clusters  are  visible  from  the  displacement  of  lattice 
atoms  near  the  cluster.  In  the  next  chapter,  formation  energies  of  small  interstitial  clusters 
and  those  of  {3 1 1 } defects  of  the  same  size  are  compared. 

5.4  Summary 

In  this  chapter,  the  formation  energy  calculations  using  MD  of  small  interstitial 
clusters  as  a function  of  cluster  size  were  discussed.  It  was  observed  that  formation 
energy  per  interstitial  decreases  with  increase  in  cluster  size  thereby  suggesting  that  there 
is  a driving  force  for  aggregation  of  interstitials  to  form  clusters.  From  the  atomic  models 
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of  these  clusters  it  was  evident  that  clusters  have  a strain  field  and  lattice  atoms  near  the 
clusters  are  displaced  from  their  equilibrium  positions. 


Ef(N)/N  (eV/interstitial) 
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Fig.  5.1  Plot  of  formation  energy  per  interstitial  (Ef(N)/N)  versus  cluster  size  (N)  for 
small  interstitial  clusters. 
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Fig.  5.2  Plot  of  binding  energy  versus  cluster  size  for  small  interstitial  clusters. 
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Fig.  5.3  Atomic  model  of  two-interstitial  cluster.  Gray  colors  indicate  strained  atoms. 
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Fig.  5.4  Atomic  model  of  three-interstitial  cluster.  Gray  colors  indicate  strained  atoms 
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Fig.  5.5  Atomic  model  of  seven-interstitial  cluster.  Gray  colors  indicate  strained  atoms. 


CHAPTER  6 

FORMATION  ENERGIES  AND  RELATIVE  STABILITY  OF  {31 1 } DEFECTS  AND 

DISLOCATION  LOOPS 

6.1  Introduction 

It  has  been  observed  that  the  duration  of  boron  transient  enhanced  diffusion 
(TED)  is  exponentially  activated  with  activation  energies  in  the  range  1.3-3. 6 eV  [Zhang 
95,  Stolk  95a],  The  activation  energy  is  related  to  the  creation,  diffusion  and  annihilation 
of  the  Si  interstitials  and  also  the  interactions  between  boron  and  boron  and  silicon 
interstitials.  Eaglesham  et  al.  [Eaglesham  94]  and  Stolk  et  al.  [Stolk  95a]  suggested  that 
boron  TED  is  a result  of  emission  of  Si  interstitials  from  extended  {311}  defects.  These 
{311}  defects  are  observed  by  Transmission  Electron  Microscopy  (TEM)  for  boron 
implantation  energies  in  the  range  less  than  1 KeV  [Downey  99]  to  energies  greater  than 
1 MeV  [Washburn  80],  It  has  also  been  observed  that  TED  can  occur  under  implant 
conditions  which  do  not  produce  {311}  defects  [Zhang  95]  which  implies  that  there  are 
submicroscopic  interstitial  clusters  such  as  boron  interstitial  clusters  (BICs)  which  release 
interstitials. 

The  {311}  defects  are  elongated  along  the  <01 1>  direction  and  their  width  is 
along  the  <23 3>  direction  and  their  habit  plane  is  the  {3 1 1 } plane.  A plan  view  weak 
beam  dark  field  TEM  image  of  the  {3 1 1 } defects  is  shown  in  fig.  6. 1 and  a high 
resolution  cross  sectional  TEM  image  is  shown  in  fig.  6.2  [Li  2000].  As  is  evident  from 
these  images  the  {311}  defects  are  rod-like  and  they  are  composed  of  chains  of 
interstitials  extending  longitudinally  along  the  <01 1>  direction,  arranged  laterally  along 
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the  <233>  direction,  and  they  lie  on  the  {3 1 1 } plane.  There  have  been  very  few 
investigations  in  the  literature  on  the  energetics  of  {3 1 1 } defects  with  due  consideration 
being  given  to  the  growth  in  width  of  the  {31 1 } defects. 

Another  aspect  of  the  {311}  defects  is  that  these  rod-like  defects  unfault  into 
dislocation  loops  upon  annealing  at  a high  temperature.  It  has  been  observed  by  TEM  that 
concomitant  with  the  dissolution  of  the  {3 1 1 } defects  some  small  loops  grow  in  size  [Wu 
77].  Similar  connection  between  the  dissolution  of  the  {3 1 1 } defects  and  the  growth  of 
loops  has  been  observed  by  other  groups  [Tan  81].  Recently,  Li  and  Jones  have  shown 
through  TEM  studies  that  sub-threshold  dislocation  loops  nucleate  from  {311}  defects  by 
unfaulting  of  the  latter  [Li  98].  This  work  involved  in  situ  annealing  plan- view 
Transmission  Electron  Microscopy  samples  in  a TEM.  It  was  observed  that  after  initial  5 
minutes  of  annealing  a dense  combination  of  both  {311}  defects  and  small  subthreshold 
dislocation  loops  was  found.  Subsequently,  when  the  sample  was  annealed  further  the 
{311}  defect  density  decreased  rapidly  and  the  loop  density  increased.  The  time 
evolution  of  500  {31 1 } defects  was  monitored.  From  this  scrutiny  it  was  observed  that 
the  unfaulting  of  a {3 1 1 } defect  was  the  source  of  every  subthreshold  loop. 

Robertson  et  al.  have  observed  the  evolution  of  both  {311}  defects  and 
dislocation  loops  for  amorphizing  implants,  in  the  end-of-range  damage  region  by  ex  situ 
TEM  [Robertson  2000],  The  post-implant  anneals  were  carried  out  in  a furnace  at  750°  C 
for  times  ranging  from  10  to  370  minutes.  It  was  seen  that  {311}  defects  were  the 
preferential  sites  for  dislocation  loop  nucleation.  However,  it  was  also  observed  that  45% 
of  the  loops  had  nucleated  by  10  min  at  750°  C. 
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In  this  chapter  a comparison  is  made  between  the  formation  energy  per  interstitial 
of  small  clusters  and  {311}  defects.  Tight-binding  calculations  performed  recently 
indicate  that  infinitely  extended  (110)  interstitial  chains  have  a formation  energy  per 
interstitial  of  1.7  eV  [Kim  97].  However,  the  formation  energy  per  interstitial  of  small 
chain  sizes  are  found  to  be  very  large  indicating  that  this  structure  is  not  suitable  for 
describing  the  early  stages  of  {31 1 } nucleation. 

Further,  tight-binding  calculations  have  shown  [Bongiorno  99]  that  there  exist 
alternative  interstitial-cluster  configurations,  which,  at  sizes  less  than  10  interstitials,  have 
formation  energy  per  interstitial  values  lower  than  that  of  the  (110)  interstitial-chain. 
Therefore,  it  can  be  concluded  from  this  data  that  small  clusters  are  not  direct  precursors 
of  the  {3 1 1 } defects  and  that  a structural  transformation  should  occur  at  some  stage  of  the 
growth  of  these  defects.  In  addition  to  this,  deep  level  transient  spectroscopy  (DLTS) 
analysis  has  shown  that  changes  occur  in  the  electrical  signatures  during  the  transition 
from  interstitial  clusters  to  {311}  defects  [Benton  97], 

More  recently,  Coffa  et  al.  [Coffa  2000]  have  used  photoluminescence  and 
transmission  electron  microscopy  (TEM)  analyses  to  study  the  transition  of  small 
interstitial  clusters  into  {311}  defects.  From  their  study  they  conclude  that  a structural 
transformation  accompanies  the  transition  from  interstitial-cluster  to  {31 1 } defects. 

Eaglesham  et  al.,  through  cross  sectional  high  resolution  electron  microscopy 
(HREM),  reported  agreement  with  the  generally  accepted  structure  of  {3 1 1 } defects 
consisting  of  chains  of  interstitials  parallel  to  the  <01 1>  direction  and  lying  on  the  {3 1 1 } 
habit  plane  with  the  chains  arranged  laterally  along  the  <23 3>  direction  [Eaglesham  94], 
Eaglesham’s  work  also  found  that,  during  Rapid  Thermal  Anneal  (RTA),  when  the  length 
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is  measured  from  plan-view  images  and  the  width  is  measured  from  HREM,  the  width 
was  found  to  be  independent  of  anneal  time  and  the  size  variations  were  increases  in 
length  only.  This  indicates  that  there  is  a saturation  of  width  after  a certain  time  of  growth 
of  the  defect. 

In  the  work  of  Cuendet  et  al.,  which  was  Monte  Carlo  simulation  work  [Cuendet 
95],  the  issue  of  elongation  and  broadening  of  the  {3 1 1 } defect  has  been  examined 
assuming  infinitely  long  stacking  fault  (SF)  ribbons  for  elongation  only  and  SF  loops  for 
both  lengthening  and  broadening.  The  defect  formation  energies  for  defect  size 
corresponding  to  60  interstitials  are  in  actuality  values  of  infinitely  long  ribbons  of  a 
particular  width.  It  was  also  noticed  in  the  work  of  Cuendet  et  al.  work  that  after  the 
ribbon  is  10  dimers  long,  the  formation  energy  per  interstitial  is  almost  constant.  It  was 
observed  that  this  formation  energy  decreases  as  the  width  widens.  However,  in  this 
work,  the  simultaneous  growth  in  length  and  width  was  not  studied  since  an  infinitely 
long  ribbon  was  assumed  for  all  the  widths  studied. 

Gencer  and  Dunham  [Gencer  97]  have  developed  a model  for  calculation  of  the 
kinetic  growth  rate  for  {3 1 1 } defects.  In  their  model  they  have  assumed  that  {311} 
defects  have  a rectangular  shape  and  that  the  {3 1 1 } defects  grow  primarily  in  their 
length.  In  their  work  it  was  assumed  that  the  {3 1 1 } defects  have  a maximum  average 
width  of  about  5nm.  It  was  also  stated  in  this  model  that  it  would  be  unreasonable  to 
assume  that  the  width  is  always  5nm,  thereby  larger  than  their  length(l)  for  1 < 5nm.  It 
was  therefore  concluded  that  the  width  (w)  may  be  expressed  as  w = min(l,5  nm).  This 
satisfies  the  condition  that  1 is  always  greater  than  or  equal  to  w and  that  w=5  nm  for 
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large  sizes  which  can  be  observed  under  high-resolution  TEM.  Here  again,  the  evolution 
of  both  length  and  width  with  increase  in  defect  size  was  ignored. 

Hobbler  and  Rafferty  [Hobler  99]  have  developed  a model  for  {31 1 } defect 
evolution  based  on  the  rate  equations  approach.  In  their  work  the  strain  energy  of  the 
{311}  defect  is  calculated  within  the  framework  of  the  theory  of  dislocations  in  isotropic 
continua  as  a special  case  of  piecewise  straight  dislocation  configurations.  The  strain 
energy  consists  of  the  self-energies  of  the  four  dislocation  segments  (two  long  edges 
along  the  sides  and  two  short  edges  along  the  ends)  and  also  the  interaction  energies  of  all 
pairs  of  these  segments.  From  this  expression,  the  variation  of  the  strain  energy  per  atom 
is  obtained  as  a function  of  defect  size  (the  number  of  interstitials)  for  2 to  6 chains  wide. 
The  results  thus  obtained  seem  to  agree  with  Cuendef  s results  excepting  at  the  smallest 
sizes  where  the  continuum  approach  is  not  valid. 

In  the  work  of  Kim  et  al.  [Kim  97]  tight-binding  calculations  were  carried  out  to 
study  the  structural  properties  and  energetics  of  extended  {311}  defects  based  on  their 
dimensions  and  interstitial  concentrations.  Their  finding  was  that  interstitial  chain 
structures  along  the  <11 0>  direction  are  energetically  more  favorable  than  isolated 
interstitials.  They  also  found  that  interstitial  chains  are  basic  building  blocks  of  the 
extended  {311}  defects.  They  suggested  that  the  extended  {311}  defects  are  formed  by 
condensation  of  the  interstitial  chains  side  by  side  in  the  <23  3>  direction.  However  their 
conclusion  was  that  the  growth  of  {3 1 1 } defects  would  occur  first  by  elongation  of 
interstitial  chains  along  the  <01 1>  direction  followed  by  an  increase  in  width  by  the 
capture  of  interstitial  chains  alongside  the  <23 3>  direction.  They  did  not  determine  the 
progress  of  defect  growth  both  longitudinally  and  laterally  at  the  same  time.  Therefore,  in 
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their  work,  they  did  not  observe  the  ultimate  saturation  in  width  which  is  observed 
experimentally.  Nevertheless  the  interstitial  configuration  developed  by  them  for 
simulation  purposes  is  found  to  be  valuable  for  carrying  out  conjugate  gradient 
calculations  as  part  of  this  thesis. 

From  the  preceding  discussion  it  is  evident  that  various  groups  have  tried  to 
examine  the  morphological  changes  in  {311}  defects  with  increase  in  defect  size.  Many 
assumptions  have  been  made  by  these  groups  in  studying  the  energetics  of  {311}  defects. 
However,  no  attempt  has  been  made,  through  atomistic  simulations,  to  elucidate  the 
simultaneous  increase  in  width  along  with  length  as  the  defect  grows.  It  has  been 
observed  experimentally  that  as  the  defect  grows  in  length  there  is  initially  a growth  in 
width  also.  However  the  width-wise  growth  is  not  indefinite  and  a saturation  is  observed 
in  the  width  after  a certain  size.  The  width  has  been  observed  to  grow  to  approximately 
60  A before  any  further  growth  stops  [Li  2000]. 

In  this  work,  for  the  first  time,  through  simulations,  the  stability  of  the  {311} 
defects  has  been  determined,  depending  on  the  length  and  width  of  the  interstitial  chains 
that  constitute  the  defect.  An  attempt  has  been  made  to  relate  the  simulation  results  to 
experimental  evidence.  An  analysis  has  also  been  done,  through  simulations,  of  the 
relative  stability  of  {311}  defects  and  perfect  dislocation  loops.  The  methodology  of 
simulations  has  been  the  Conjugate  Gradient  method  using  the  Stillinger-Weber 
interatomic  potential. 

6.2  Computational  Details 

As  mentioned  in  previous  chapters,  the  Stillinger  -Weber  (SW)  interatomic 
potential  is,  by  far,  the  most  widely  used  interatomic  potential  for  silicon  and  therefore  is 
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the  potential  of  choice  in  this  work.  The  Conjugate  Gradient  technique  is  used  in  this 
investigation.  This  technique  is  also  called  Molecular  Statics  and  unlike  in  Molecular 
Dynamics(MD)  the  trajectories  of  atoms  are  not  monitored.  The  energy  minimization  is 
carried  out  at  0 K.  Essentially  the  system  is  relaxed  at  OK  by  minimizing  the  potential 
energy  function.  The  computational  time  is  vastly  less  than  the  MD  technique  and 
therefore  large  systems  can  be  used  in  this  technique. 

6.2.1  The  Mathematics  Behind  the  Conjugate  Gradient  Method 

The  case  of  calculating,  at  a given  N-dimensional  point,  not  just  the  value  of  a 
function  f(P)  but  also  the  gradient  Vf(P)  has  to  be  considered  first.  If  the  gradient  is  not 
used  N2  separate  line  minimizations  have  to  be  made  in  calculating  the  function 
minimum,  each  requiring  a few  function  evaluations.  Now  each  evaluation  of  the  gradient 
will  bring  N new  components  of  information.  If  used  judiciously,  we  need  to  make  only 
of  order  N line  minimizations. 

A factor  of  N improvement  in  computational  speed  in  not  necessarily  implied.  As 
a rough  estimate,  the  calculation  of  each  component  of  the  gradient  takes  about  as  long  as 
evaluating  the  function  itself.  In  that  case  there  will  be  of  order  N2  equivalent  function 
evaluations  both  with  and  without  gradient  information.  Even  if  the  advantage  is  not  of 
order  N,  it  is  quite  substantial.  This  is  due  to  the  following  reasons:  (a)  Each  computed 
component  of  the  gradient  information  will  not  just  save  one  function  evaluation  but  a 
number  of  them,  equivalent  typically  to,  a whole  line  minimization,  (b)  There  is  usually  a 
high  degree  of  redundancy  in  the  formulas  for  various  components  of  a function’s 
gradient:  if  this  is  the  case,  especially,  when  there  is  also  redundancy  with  the  calculation 


78 


of  the  function,  then  the  calculation  of  the  gradient  may  cost  significantly  less  than  N 
function  evaluations. 

However,  it  is  erroneous  to  assume  that  any  reasonable  way  of  incorporating 
gradient  information  should  be  about  as  good  as  any  other.  The  steepest  descent  method 
is  an  example  of  this  not  very  good  algorithm.  The  steepest  descent  method  starts  at  a 
point,  say  P0,  and  as  many  times  as  necessary  moves  from  point  Pj  to  the  point  P,+1  by 
minimizing  along  the  line  from  Pi  in  the  direction  of  the  local  downhill  gradient  -Vf(P,). 
The  method  will  perform  many  small  steps  in  going  down  a long,  narrow  valley,  even  if 
the  valley  is  a perfect  quadratic  form.  It  may  be  expected  that  the  first  step  would  take  us 
to  the  local  minimum  and  the  subsequent  step  would  be  down  the  long  axis  of  the  valley. 
Since  the  new  gradient  at  the  local  minimum  point  of  any  line  minimization  is 
perpendicular  to  the  direction  just  traversed,  with  the  steepest  descent  method,  a right 
angle  turn  must  be  made  which  does  not,  in  general,  take  us  to  the  minimum  of  the 
function. 

What  is  really  required  is  to  travel  not  down  the  new  gradient,  but  rather  in  a 
direction  that  is  somehow  conjugate  to  the  old  gradient,  and  as  far  as  possible,  to  all  the 
previously  traveled  directions.  Methods  that  achieve  this  construction  are  called 
conjugate  gradient  methods.  If  a function  f(x)  can  be  approximated  as 

f (x)  = c - b.x  + (1/2)  x.A.x  (6.1) 

where  the  matrix  A whose  components  are  the  second  partial  derivative  matrix  of  the 
function  is  called  the  Hessian  matrix  of  the  function  at  some  particular  point  P. 


Using  this  approximation  of  equation  (6.1),  the  gradient  of  the  f is  calculated  as 

Vf=  A.x-band6(Vf)  = A.(8x)  (6.2) 
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This  implies  that  the  function  will  be  at  an  extremum  at  a value  of  x obtained  by  solving 
A.x  = b.  Suppose  that  we  have  moved  along  some  direction  u to  a minimum  and  now 
decide  to  move  along  some  new  direction  v.  Then,  the  condition  that  motion  along  v not 
spoil  our  minimization  along  u is  just  that  the  change  in  the  gradient  be  perpendicular  to 
u.  This  is  the  basis  for  the  conjugate  gradient  technique.  From  equation  (6.2)  this  is  given 
as  follows: 

u.8(  Vf)  = u.A.v  = 0 (6.3) 

When  equation  (6.3)  holds  for  two  vectors  u and  v they  are  said  to  be  conjugate.  When 
this  equation  holds  pair-wise  for  all  members  of  a set  of  vectors,  they  are  said  to  be  a 
conjugate  set.  If  we  carry  out  successive  line  minimization  of  a function  along  a 
conjugate  set  of  directions,  then  we  need  not  redo  any  of  those  directions.  If  a direction 
set  of  N linearly  independent  mutually  conjugate  directions  is  created,  then,  one  pass  of  N 
line  minimizations  will  put  it  exactly  at  the  minimum  of  a quadratic  form  such  as  that 
given  by  equation  (6.1). 

6.2.2  Interstitial  Configuration  for  Calculations 

In  this  work  a 27000  atom  computational  cell  was  used.  This  large  cell  size 
ensures  that  the  lattice  displacement  field  created  by  the  relaxed  extended  defects  is 
contained  within  the  cell.  For  {311}  defects,  the  interstitial  configuration  used  here  is 
similar  to  that  used  by  Kim  et  al  [Kim  97],  This  model  is  illustrated  in  fig.  6.3.  The 
interstitial  chains  constituting  the  {31 1 } defect  are  extended  along  the  <01 1>  direction, 
laterally  along  the  <233>  direction  and  the  chains  lie  on  the  {311}  plane.  The  adjacent 
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interstitial  chains  are  separated  from  each  other  by  a distance  of  a — where  a is  the 

2V2 

lattice  parameter.  This  configuration  is  then  relaxed  using  the  conjugate  gradient  method 
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with  the  SW  potential.  As  a consequence  of  this  relaxation,  the  formation  energy  of  the 
{311}  defects  is  obtained  as  a function  of  the  number  of  interstitials  in  the  defect. 
Calculations  are  done  for  different  number  of  chains  up  to  a configuration  five  chains 
wide.  Similarly,  the  interstitial  configuration  of  the  perfect  dislocation  loop  with  Burgers 
vector  a/2<l  10>  is  created  and,  subsequently,  the  loop  configuration  is  relaxed  using  the 
conjugate  gradient  technique  to  yield  formation  energies  for  varying  loop  size. 

6.3  Results  and  Discussion 

First,  fig.  6.4  is  a plot  of  the  formation  energies  of  small  clusters  and  {31 1 }s  of 
the  same  size  and  it  shows  the  difference  between  the  formation  energies  of  both  types  of 
defects.  The  cluster  formation  energies  were  obtained  in  the  previous  chapter  from 
molecular  dynamics  calculations  and  the  {31 1 } formation  energies  are  obtained  in  this 
chapter  from  conjugate  gradient  calculations.  Thus,  these  results  are  similar  to  those 
obtained  from  tight-binding  results  mentioned  in  a previous  section  of  this  chapter 
[Bongiorno  99],  It  may  be  concluded  that  the  small  interstitial  clusters  are  not  direct 
precursors  of  {31 1 }s  and  that  clusters  undergo  a structural  transformation  during  their 
growth  before  they  form  {311}  defects. 

The  total  formation  energies  for  the  {31 1 } defects  obtained  by  the  Conjugate 
Gradient  method  is  plotted  as  a function  of  defect  size  in  fig.  6.5.  It  is  evident  from  this 
plot  that  as  the  defect  size  increases,  the  least  energy  configuration  has  progressively 
increasing  widths.  This  trend  is  shown  for  the  first  four  chains  only.  The  five-chain  wide 
defect  has  a higher  formation  energy  compared  to  the  four-chain  defect  for  the  defect 
sizes  accessible  from  this  simulation.  However,  a closer  inspection  of  the  data  seems  to 
suggest  that  the  five-chain  defect  will  become  more  stable  than  the  four-chain  wide 
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defect  at  a larger  defect  size.  Moreover,  it  is  also  indicated  from  this  plot  that  the  total 
formation  energy  varies  almost  linearly  with  defect  size. 

Figure  6.6  shows  the  formation  energy  per  interstitial  of  the  {311}  defect  as  a 
function  of  the  defect  size  for  different  number  of  chains  that  constitute  the  defect.  This 
plot  indicates  that,  for  the  defect  sizes  studied,  the  formation  energy  per  interstitial 
decreases  as  the  defect  size  increases.  This  suggests  that  there  is  a driving  force  for  defect 
growth.  Moreover,  as  in  fig.  6.5,  the  trend  shown  is  that  the  most  stable  defect  width 
(number  of  chains)  increases  as  the  defect  size  increases.  However,  it  is  also  seen  that  this 
increase  in  defect  width  does  not  continue  indefinitely  since  the  five-chain  wide  defect  is 
less  stable  than  the  four-chain  wide  defect.  Therefore,  there  seems  to  be  a saturation  in 
the  defect  width  with  increasing  defect  size. 

Figure  6.7  illustrates  this  point  further.  It  is  a plot  of  the  formation  energy  per 
interstitial  as  a function  of  width  (number  of  chains)  for  a defect  size  of  30  interstitials.  It 
is  seen  in  this  plot  that  for  this  defect  size  of  30  interstitials  the  most  stable  defect  has  a 
width  of  3 chains.  Also,  the  slopes  of  this  curve  on  either  side  of  this  minimum  point  are 
different.  On  the  right  of  this  minimum  the  slope  is  much  smaller  than  that  on  the  left. 
This  implies  that  the  width  is  beginning  to  saturate  as  the  defect  size  increases.  This 
agrees  well  with  our  earlier  conclusion  that  the  width  does  not  increase  indefinitely  with 
increase  in  defect  size. 

It  has  been  observed  experimentally  [Li  2000]  that  the  average  width  of  the  {3 1 1 } 
defect  initially  increases  with  average  length  and  then  saturates  at  a value  of  around  60  A. 
This  observation  was  made  using  cross-sectional  HREM.  This  trend  is  in  accordance  with 
the  energetics  calculation  done  in  this  work.  Fig.  6.8  is  a plot  of  {31 1 } width  versus 
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{311}  length  observed  experimentally  [Li  2000]  as  well  as  a length  to  width  correlation 
obtained  from  these  simulations.  In  order  to  obtain  the  average  length  for  each  width 
from  the  simulation  data,  a certain  width  size  is  chosen  and  the  average  defect  size  during 
which  this  width  is  most  stable  is  measured  from  the  energetics  plot  (fig. 6.6).  Upon 
dividing  the  average  size  by  the  number  of  atoms  in  width  the  average  number  of  atoms 
in  length  is  obtained.  From  this  the  average  length  of  the  defect  is  estimated  using  the 
model  of  Kim  et  al.  [Kim  97]  which  is  discussed  in  section  6.2.2.  This  data  is  computed 
for  three  chain-wide  and  four  chain-wide  defects  and  shows  a direct  correlation  of  length 
versus  width.  In  the  plot  of  experimental  data  the  saturation  of  the  width  is  evident.  In 
order  to  confirm  the  validity  of  the  trends  seen  from  the  atomistic  simulations  presented 
in  this  thesis,  a calculation  was  done  based  on  dislocation  theory  of  the  width  of  the 
stacking  fault  defect.  This  calculation  was  done  as  follows: 

The  maximum  separation  between  two  partial  dislocations  is  dictated  by  the 
stacking  fault  energy  of  the  stacking  fault  separating  the  two  dislocations.  The  maximum 
width  of  the  stacking  fault  is  given  as  w in  the  following  equation: 

pb2/(27tw)  = y (6.4) 

where  p is  the  shear  modulus  of  silicon  (p  = 6.5x1 010  Pa),  b is  the  length  of  the 
Burgers  vector  (b  = 0.1071  nm),  and  y is  the  stacking  fault  energy  in  silicon  (y=  60 
mJ/m2) . From  this  calculation  the  width  of  the  stacking  fault,  w = 19.8  A.  This  is  in 
agreement  with  the  width  obtained  using  atomistic  simulations  which  has  a value  of  20 

A. 
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The  driving  force  for  the  formation  of  defects  such  as  the  {311}  defects  and 
dislocation  loops  is  the  minimization  of  the  free  energy  of  the  system.  The  free  energy  of 
a defect  of  size  n is  given  by 

AG„  = -nkT  ln(C,/C*)  + AGnf  (6.5) 

Here,  Ci/C*  is  the  initial  supersaturation  of  interstitials  in  the  lattice  and  AG,/  is 
the  formation  energy  of  defect  of  size  n.  C*,  the  equilibrium  concentration  of  interstitials, 
is  calculated  as  described  in  chapter  4 using  the  interstitial  formation  energy  (for  the 
Stillinger-Weber  interatomic  potential)  and  formation  entropy.  These  calculations  are 
done  at  750  °C.  Defect  formation  energy  is  obtained  from  conjugate  gradient  calculations 
as  described  in  previous  sections  of  this  chapter.  Formation  energy  data  is  available  for 
both  {311}  defects  and  perfect  dislocation  loops  through  conjugate  gradient  calculations. 
Plots  of  the  free  energy,  AG„,  as  a function  of  defect  size  n yield  nucleation  energy  barrier 
values  and  critical  nucleus  size  values  for  the  nucleation  process. 

Figs.  6.10-6.15  show  plots  of  free  energy  versus  defect  size  for  {31 1 }s  and 
perfect  loops  at  various  doses.  The  interstitial  concentration  (Ci)  is  obtained  from 
supersaturation  values  (Ci/C*).  The  value  of  C*  is  estimated  by  the  technique  used  in 
chapter  4.  The  formula  for  estimating  C*  is  as  follows: 

C*  = 5xl022  exp[-(Ef-  TSf)/kBT]  /cm3  (6.6) 

where  Ef  is  the  formation  energy  of  the  interstitial  which  is  obtained  from  molecular 
dynamics  simulations  as  described  in  chapter  4 and  T is  the  absolute  temperature.  The 
formation  entropy  (Sf)  is  estimated  by  the  approach  adopted  by  Tang  et  al.  [Tang  97]  as 
described  in  chapter  4.  The  value  of  Ci  is  obtained  by  multiplying  the  supersaturation 
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(Ci/C*)  by  C*.  The  value  of  the  implant  dose  0 is  estimated  by  using  the  following 
expression  which  relates  dose  to  peak  concentration,  in  this  case  for  a self  implant 
Ci  = (0.4  x 0)/ ARp  (6.7) 

where  ARpis  the  ion  projected  straggle  for  silicon  ions  which  is  almost  equal  to  that  of 
phosphorus  ions  which  in  turn  is  available  in  textbooks  on  silicon  processing.  From  these 
plots  it  is  evident  that,  around  a dose  of  1 .2x1 012/cm2,  free  energy  curves  for  the  {3 1 1 } 
defects  begin  to  roll  over  indicating  that  the  {311}  defects  are  beginning  to  nucleate.  Also 
from  these  plots,  considering  the  nucleation  barrier,  at  a dose  of  around  8xl0l3/cm2,  there 
is  a cross-over  from  {311}  defects  to  loops  as  far  as  stability  is  concerned.  This  is 
illustrated  in  fig.  6.16  where  the  barrier  height  is  plotted  as  a function  of  implant  dose. 
However,  it  should  be  noted  that  the  critical  nucleus  size  of  {3 1 1 } defects  becomes  lower 
than  that  of  the  loops  as  the  dose  increases  beyond  lxlOl4/cm2and  the  barrier  heights  are 
not  vastly  different  when  the  dose  increases  beyond  lxl014/cm2.  It  may  therefore  be 
inferred  that  the  nucleation  of  {31 1 } defects  is  possible  with  smaller  agglomeration  of 
interstitials  at  doses  larger  than  lxlOl4/cm2  since  the  critical  nucleus  size  of  these  defects 
is  smaller  than  that  of  the  loops. 

For  example,  fig.6.14  shows  a typical  case  of  free  energy  of  nucleation 
calculations  for  a dose  greater  than  lxl014/cm2.  It  is  evident  in  this  plot  that  the 
nucleation  barriers  are  not  significantly  different  although  the  critical  nucleus  size  for  the 
{311}  defect  is  6 atoms  large  whereas  for  the  loop  it  is  10  atoms  large.  As  the  interstitial 
agglomerates  form,  early  in  the  agglomeration,  it  becomes  easier  for  the  interstitials  to 
nucleate  the  {3 1 1 } defect  since  the  {3 1 1 } defect  has  a smaller  critical  nucleus  size  as 
compared  to  the  loop.  This  would  explain  the  experimental  evidence  that  at  higher  doses 
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{31 1 }s  nucleate  first  and  then  these  defects  unfault  to  form  loops.  Also,  experimentally, 
the  critical  dose  for  {31 1 } nucleation  is  7xl012/cm2  and  loops  become  more  stable  than 
{3 1 1 }s  at  a dose  around  2x1 0l4/cm2  [Jones  96],  Therefore,  we  find  that  the  theoretical 
calculations  carried  out  in  this  work  yield  trends  which  conform  to  experimental 
evidence.  Figs.  6.10-6.16  support  this  conclusion. 

According  to  Li  et  al.  [Li  98],  the  {31 1 }s  are  the  source  of  the  loops  and  this 
experiment  was  done  at  sub-amorphizing  implant  doses.  Fig.  6.9  is  a plot  of  the 
formation  energy  per  interstitial  as  a function  of  defect  size  for  {31 1 } defects  and  perfect 
loops.  From  this  plot  it  is  evident  that  for  both  {3 1 1 }s  and  loops  the  formation  energy  per 
interstitial  decreases  continuously  with  defect  size.  This  implies  that  there  is  a driving 
force  for  increase  in  defect  size.  The  implications  of  the  trends  in  defect  nucleation  can  be 
validated  with  an  inspection  of  the  experimental  data  obtained  by  Robertson  et  al. 
[Robertson  2000]  for  amorphizing  implants. 

Fig.  6.17  is  a plot  of  percentage  loop  nucleated  as  a function  of  anneal  time  for  the 
experiment  [Robertson  2000]  described  in  the  first  part  of  this  chapter.  In  that  work,  it 
was  noticed  that  45%  of  the  loops  were  nucleated  in  the  first  10  minutes  of  anneal.  Out  of 
these,  the  nucleation  of  25%  of  the  loops  could  not  be  explained  as  having  occurred  due 
to  unfaulting  of  {3 1 1 } defects.  This  number  of  unexplained  loops  was  determined  by  the 
fact  that  they  did  not  nucleate  at  the  original  location  of  {3 1 1 }s.  By  overlapping  of 
sequential  TEM  micrographs  before  and  after  unfaulting,  it  was  found  that  75%  of  the 
nucleated  loops  were  located  at  the  same  spots  as  that  of  the  original  {3 1 1 }s  and  the 
remaining  25%  could  not  be  found  to  originate  at  {311}  sites. 
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The  free  energy  calculations  carried  out  in  this  chapter  show  that  the  {3 1 1 } 
defects  nucleate  first  during  an  anneal.  The  unfaulting  reaction  to  convert  {311}  defects 
into  perfect  dislocation  loops  is  given  as  follows  [Li  98]: 

^-[119+  ^[19  2 9]  = |[  101]  (6.8) 

As  an  example,  consider  a {31 1 } defect  of  length  1 100 A and  width  45 A as  observed 
experimentally  by  Jinghong  et  al.  [Li  2000],  The  total  number  of  atoms  contained  in  this 
defect  is  estimated  to  be  1710  using  the  {311}  interstitial  configuration  model  developed 
by  Kim  et  al.  [Kim  97]  which  is  described  in  section  6.2.2.  The  total  energy  of  this  {311} 
defect  is  calculated  by  estimating  Gb2  x total  length  of  the  {31 1 } + (length)x(width)x(y) 
where  G is  the  shear  modulus  of  silicon  (6.45x10 10  Pa),  b is  the  Burgers  vector  of  the 
{311}  defect  (a/21  [1 16])  and  y is  the  stacking  fault  energy  of  silicon  (60  mJ/m2).  This 
total  energy  is  estimated  from  the  calculation  to  be  2532  eV.  Correspondingly,  the  radius 
(R)  of  the  final  loop  is  estimated  to  be  60A  using  the  silicon  atom  density  on  the  { 1 1 1 } 
plane  (1.5xlOl5/cm2)  with  the  total  number  of  atoms  in  the  loop  being  1710  using  the 
expression  7tR2x(atom  density)  = Total  number  of  atoms  . The  total  energy  of  the  loop  is 
calculated  to  be  27tRGbi2  where  b|  is  the  Burgers  vector  of  the  loop  (a/2[  1 01  ]).  This  total 
energy  value  for  the  loop  is  calculated  to  be  2241  eV.  The  energy  barrier  for  the 
unfaulting  reaction  is  estimated  by  computing  the  energetics  of  the  transform  partial 
which  sweeps  through  the  stacking  fault  during  the  unfaulting  reaction.  This  has  a value 
of  Gbp2  x(length)  = 3300eV,  where  bp  is  the  Burgers  vector  of  the  transform  a/42[19  2 9] 
. The  {311}  defect  unfaults  into  a loop  upon  annealing  the  silicon  sample  at  high 
temperatures  such  as  750  °C  which  enables  the  {31 1 } defect  to  surmount  the  energy 
barrier  of  unfaulting.  The  25%  of  unexplained  loops  in  the  work  of  Robertson  et  al. 
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[Robertson  2000]  could  originate  in  submicroscopic  interstitial  defects  which  are 
primitive  {311}  defects  and  this  is  supported  by  nucleation  theory  results  explained 
above. 


6.4  Summary 

It  was  observed  that,  for  cluster  size  less  than  1 0 interstitials,  formation  energy 
per  interstitial  values  were  lower  than  that  of  {3 1 1 } defects  in  the  same  size  regime 
implying  that  small  interstitial  clusters  are  not  direct  precursors  of  {31 1 } defects.  It  is 
inferred  that  the  small  clusters  have  to  undergo  a structural  transformation  before 
nucleating  {311}  defects.  Atomistic  simulations  have  been  done  using  the  conjugate 
gradient  method  with  the  Stillinger-Weber  potential  to  obtain  the  formation  energy  of 
{311}  defects  for  defect  widths  up  to  5 chains. 

It  is  observed  that  the  width  of  the  most  stable  defect  increases  initially  at  first 
with  defect  size  and  then  saturates.  This  trend  is  in  agreement  with  experiment.  The 
formation  energy  of  perfect  loops  was  also  obtained  by  conjugate  gradient  simulations. 
From  free  energy  calculations  it  is  found  that  {311}  defects  nucleate  earlier  in  the 
annealing  phase  and  unfault  into  loops  as  the  anneal  proceeds.  This  explains  the 
nucleation  of  loops  in  the  end  of  range  region  for  amorphized  and  annealed  samples. 
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Fig.  6.1  Plan  view  weak  beam  dark  field  image  of  {31 1}  defects. 
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Fig.  6.2  Cross-section  HRTEM  image  of  {3 1 1 } defect. 
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Fig.  6.3  Configuration  of  {31 1 } defects.  The  dimensions  of  the  unit  cell  are  Lx0  = a/V2, 
Ly0  = a Vl  1/(2V2)  and  Lz0  = aVl  1,  where  a is  the  Si  lattice  constant.  The  size  of  the 
computational  cell  is  Lx  = nx  x Lx0,  Ly  = ny  x Ly0  and  Lz  = 2 x Lz0. 
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Fig.  6.4  Plots  of  formation  energy  per  interstitial  (Ef(N)/N)  as  a function  of  defect  size(N) 
for  small  clusters  and  {311}  defects  of  corresponding  size. 
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Fig.  6.5  Total  formation  energy  (Ef(N))  of  {3 1 1 } defect  as  a function  of  defect  size  for 
different  widths 
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Fig.  6.6  Formation  energy  per  interstitial  (Ef(N)/N)  for  {311}  defect  as  a function  of 
defect  size  for  different  widths. 
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Fig.  6.7  Formation  energy  per  interstitial  (Ef(N)/N)  as  a function  of  width  (number  of 
chains)  for  a defect  size  of  30  interstitials. 
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• Width  of  {311}  vs  length  of  {31 1}(Experimental) 
— * — Width  of  {311}  vs  length  of  {311}  (Simulations) 
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Fig.  6.8  Width  as  a function  of  average  length  of  {3 1 1 } defects 
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Fig.  6.9  Formation  energy  per  interstitial  (Ef(N)/N)  for  {3 1 1 }s  (for  different  widths)  and 
perfect  dislocation  loop 
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Fig.  6. 10  Free  energy  of  nucleation  as  a function  of  defect  size  for  the  {3 1 1 } defect  at  a 
dose  of  8x1 01  '/cm2 
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Fig.  6. 1 1 Free  energy  of  nucleation  as  a function  of  defect  size  for  the  {3 1 1 } defect  at  a 
dose  of  1.2xlOl2/cm2 
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Fig.  6. 12  Free  energy  of  nucleation  as  a function  of  defect  size  for  the  {3 1 1 } defect  and 
loop  at  a dose  of  4.7x1 012/cm2 
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Fig.  6. 1 3 Free  energy  of  nucleation  as  a function  of  defect  size  for  the  {3 1 1 } defect  and 
loop  at  a dose  of  1.2x1 0l3/cm2 
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Fig.  6.14  Free  energy  of  nucleation  as  a function  of  defect  size  for  the  {311}  defect  and 
loop  at  a dose  of  1.2x1 014/cm2 
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Fig.  6.15  Free  energy  of  nucleation  as  a function  of  defect  size  for  the  {311}  defect  and 
loop  at  a dose  of  2.4xl014/cm2 
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Fig.  6.16  Nucleation  barrier  height  as  a function  of  implant  dose  for  the  {311}  defect  and 
loop 
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Fig.  6.17  Percentage  loops  nucleated  as  a function  of  annealing  time 


CHAPTER  7 
CONCLUSIONS 

Ion  implantation  will  continue  to  be  the  most  widely  used  technique  of 
introducing  dopants  into  silicon  substrates  in  integrated  circuit  manufacturing.  Ion 
implantation  produces  damage  and  upon  post-implantation  annealing  the  excess 
interstitials  agglomerate  to  form  submicroscopic  clusters  and  extended  defects.  It  is  now 
well  established  that  point  defects,  submicroscopic  clusters  and  extended  defects  play  an 
important  role  in  dopant  diffusion.  This  is  particularly  important  while  studying  the 
phenomenon  of  transient  enhanced  diffusion  which  is  a limiting  factor  in  the  formation  of 
shallow  junctions. 

Since  point  defects  have  to  diffuse  in  the  lattice  before  agglomerating  to  form 
clusters,  it  is  important  to  analyze  the  point  defect  diffusivities  and  self  diffusion.  While 
experimental  methods  cannot  yield  information  about  diffusivities  of  individual  point 
defects  such  as  interstitials,  vacancies  and  also  the  di-vacancy,  theoretical  calculations 
have  been  found  to  be  very  useful  in  this  task.  Particularly,  molecular  dynamics 
simulations  can  be  used  to  calculate  the  point  defect  diffusivities  since  it  is  possible  using 
this  technique  to  compute  mean  squared  displacements  over  all  atoms  in  the  presence  of 
defects.  The  interatomic  potential  used  in  molecular  dynamics  is  crucial  as  far  as  the 
validity  of  the  results  are  concerned.  Among  the  over  30  interatomic  potentials  for 
silicon  available  in  the  literature,  the  Stillinger-Weber  potential  is  by  far  the  most  widely 
used  interatomic  potential. 


105 


106 


Recently,  a new  interatomic  potential  called  the  environment-dependent 
interatomic  potential  (EDIP)  was  developed,  whose  authors  claimed  that  it  was  superior 
to  the  conventional  potentials.  In  chapter  four  of  this  dissertation  a comparative  study  of 
diffusivities  of  point  defects  and  the  di-interstitial  was  done  using  the  Stillinger-Weber, 
Tersoff  and  EDIP  potentials.  From  the  migration  energies  thus  obtained  it  was  found  that 
the  EDIP  potential  lacks  validity  as  it  yields  a very  low  value  for  migration  energy  of  the 
interstitial  (0.334  eV)  as  compared  to  the  value  given  by  the  other  two  potentials  as  well 
as  results  from  LDA  available  in  the  literature  (0.9-1 .0  eV).  Moreover,  using  formation 
and  migration  energies  the  self-diffusion  coefficient  was  calculated  for  the  interstitial  and 
vacancy.  It  was  found  that  the  crossover  temperature  for  interstitials  and  vacancy  self- 
diffusion  coefficients  was  717  °C  whereas  the  experimental  value  is  1050  °C  which 
further  proves  that  the  EDIP  potential  is  flawed. 

In  chapter  5 of  this  dissertation,  the  results  were  presented  of  molecular  dynamics 
simulations  with  the  Stillinger-Weber  potential  used  to  estimate  the  formation  energy  of 
small  interstitial  (<10  interstitials)  clusters.  It  was  found  that  the  formation  energy  per 
interstitial  decreases  with  increasing  cluster  size  initially  thereby  indicating  a driving 
force  for  the  aggregation  of  interstitials  into  clusters.  It  was  also  found  from  atomic 
models  that  the  clusters  have  a strain  field  in  their  vicinity.  Moreover,  in  chapter  6,  a 
comparison  between  the  formation  energy  per  interstitial  of  the  submicroscopic  clusters 
and  {311}  defects  in  the  same  size  range  indicates  that  the  clusters  are  more  stable  than 
the  {311}  defects.  This  implies  that  small  interstitial  clusters  are  not  direct  precursors  of 
{311}  defects  and  that  the  small  clusters  have  to  undergo  a structural  transformation 
before  nucleating  {311}  defects. 
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In  chapter  6,  a study  is  presented  of  Conjugate  Gradient  calculations  with  the 
Stillinger- Weber  potential  to  obtain  formation  energies  of  {3 1 1 } defects  for  defect  widths 
up  to  five  chains.  It  was  observed  that  the  {311}  width  increases  as  the  defect  size 
increases  initially  and  then  saturates.  This  is  in  agreement  with  experimental 
observations.  This  contradicts  some  of  the  studies  in  the  literature  where  linear  growth 
takes  precedence  over  increase  in  width.  The  formation  energy  of  perfect  dislocation 
loops  was  also  calculated  as  a function  of  defect  size.  From  the  formation  energies  the 
free  energy  of  nucleation  was  calculated  for  {311}  defects  and  loops  for  various  implant 
doses.  This  provides  an  explanation  for  loop  nucleation  through  unfaulting  of  {311} 
defects  which  form  earlier  in  an  anneal. 
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